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Abstract for Part I 


Earthquake energy calculations are generally made through an empirical applica- 
tion of the familiar Gutenberg-Richter energy-magnitude relationships. The precise 
physical significance of these relationships is somewhat uncertain. We make use 
here of the recent increases in knowledge about the earthquake source to place 
energy measurements on a sounder physical basis. For a simple trapezoidal far-field 
displacement source-time function with a ratio x of rise time to total duration Tq, the 


seismic energy E is proportional to 


U § 


where Mo is seismic moment. As 


x(i -z) e r 0 3 

long as x is greater than 0.1 or so, the effect of rise time is not important, “ he 
dynamic energies thus calculated for shallow events are in reasonable agreement 
with the estimate E w (5X1 0 5) A/ 0 based on elastostatic considerations. Deep 
events, despite their possibly different seismological character, yield dynamic ener- 
gies which are compatible with a static prediction similar to that for shallow events. 
Studies of strong-motion velocity traces obtained near the sources of the 1971 San 
Fernando, 1966 Parkfieid, and 1979 Imperial Valley earthquakes suggest that even 
In the distance range of 1-6 km., most of the radiated energy is below 1-2 Hz. in fre- 
quency. Far field energy determinations using long period WWSSN instruments are 
probably not in gross error despite their bandlimited nature. The strong motion record 
for the Intermediate depth Bucharest earthquake of 1977 also suggests little telese- 


Ismic energy outside the pass-band of a long period WWSSN instrument. 
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Abstract for Part II 


The pattern of seismicity as a function of depth in the world, and the orientation 
of stress axes of deep and intermediate earthquakes, are explained using viscous 
fluid models of subducting slabs, with a barrier In the mantle at 670 km. 670 km is 
the depth of a seismic discontinuity, and also the depth be!ow which earthquakes do 
not occur. The barrier in the models can be a viscosity increase of an order of magni- 
tude or more, or a chemical c.iscontlnuity where vertical velocity Is zero. Log N 
versus depth, where N is the number of earthquakes, shows (1) a linear decrease to 
about 260-300 km depth, (?) a minimum near that depth, and (3) an increase 
thereafter. Stress magnitude in a subducting slab versus depth, for a wide variety of 
models, shows the same pattern. Since there is some experimental evidence that N 
is proportional to e ka , where it Is c constant and a is the stress magnitude, the 
agreement Is encouraging. In addition, the models predict down-dip compression in 
tho slab at depths below 400 km. This has been observed in earlier studies of earth- 
quake stress axes, and we have confirmed it via a survey of events occurring since 
10 77 which have been analyzed by moment tensor Inversion. At intermediate depths, 
the models predict an approximate but not precise state of down-dip tension when 
the slab is dipping. Tne observations do not show an unambiguous state of down-dip 
tension at intermediate depths, but in the majority of regions the state of stress is 
decidedly closer to down-dip tension than It is to down-dip compression. Chemical 
discontinuities above 670 km, or phase transitions with an elevation of the boundary 


in the slab, predict, when incorporated into the models, stress peaks which are not 
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mlrrored In the profile of seismicity versus depth. Models with an astheriosphere and 
mesosphere of appropriate viscosity can not only explain the state of stress 
observed In double Benioff zones, but also yield stress magnitude profiles consistent 
with observed seismicity. Models where a nonlinear rheology Is used are qualitatively 
consistent with the linear models. 



iXr 


Table of Contents 


Part I: Tbs Energy Release in Earthquakes 

1 . Introduction 

2. Dynamic Energy from Source-Time Function — _ 

3. A Simplified Procedure for Modelng Deep Focus Events 

4. Comparison with Static Energy Estimates 

6- Near Source Energy Studies and the Question 


of Frequency Content 

0. Energy and Magnitude 

7 . Conclusions 

References 



Part II: Subduction Zone Seismicity and Stress in Slabs 


3 

4 
13 
19 

29 

38 

44 

45 


1. Introduction — 

2. Observations •»•••»•• .•*•••*•••< 

3. Calculations of Stress In Subducting Slabs — 

4. Summary and Conclusions 

References 


62 

53 

100 

199 

202 


Appendix A: Mschanisms of Large intermediate and Deep 


Earthquakes, 1978-1 981 .. - 213 

Appendix B: A Brief Outline of Directional Statistics - — . 229 

Appendix C: Plots of Seismicity versus Depth for the 

World's Subduction Zones 234 



-x- 


Appendlx D: Spatial S trass Plots and Velocity Plots 


not Included in the Text 


•••••••••••••• •»••• 




206 



- 1 - 
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Abs tract 


Earthquake energy calculations are generally made through an empirical applica- 
tion of the familiar <autenberg-Richt«r energy-mr jnitude relationships. The precise 
physical significance of these relations nips is somewhat uncertain. Wa make use 
here of the recent Increases In knowledge about the earthquake source to place 
energy measurements on a sounder physical basis. For a simple trapezoidal far-held 
displacement source-time function with e ratio x of rise time to total duration i v'e 


seismic < energy 


is proportional to 


x(1 -x) 8 


Uti 

Tg' 


when* 


is seismic moment. As 


long as x la greater than 0.1 or so, the effect of rise time is not important. The 
dynamic energies thus calculated for shallow events are in reasonable a&eement 
with the estimate E « (6X1G" 8, a/ 0 based on elastostakic considerations. Deep 
events, despite their possibly different seismologies! character, yield dynamic ener- 
gies which are compatible with a static prediction similar to that for shallow 6v *nts. 
Studies of strong-motion velocity traces obtained near the sources of the 1971 San 
Fernando, 1906 Parkfield, and 1979 Imperial Valley earthquakes suggest that ever 
In the distance range of 1-6 km., ^>ost of the radiated energy is below 1-2 Hz. in fre- 
quency. Far field energy de terminations using long period WWSSN instruments are 
probably not In gross error despite their bandlUnited nature. The strong motion record 
for the Intermediate depth Bucharest earthquake of 197? also suggests little teh'se- 
Ismic energy outside the pass-band of a long period V/WSSN instrument. 
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1. Introduction 


The energy released in earthquakes can be estimated in a number of ways (for a 
comprehensive review see Bath, 1 966). We may divide the energy estimates from 
the variety of methods available into two broad classes: the static estimates and the 
dynamic estimates. Static estimates can be obtained from static values of moment 
and stress drop; dynamic estimates, on the other hand, are obtained from seismo- 
grams. 

We review static estimates of energy In Section 4. We discuss there that with 
some simple assumptions, a static estimate of energy can be obtained from the for- 
mula E = (5X1O" s >A/ 0 (Knopoff, 1968; Kanamori, 1977). 

We may subdivide dynamic estimates of energy from body waves into two 
groups. One procedure involves the direct integration of an observed waveform at a 
particular station; another Involves integration of an inferred displacement source- 
time function. 

The familiar energy-magnitude relationships of Gutenberg and Richter (1942, 
1966a, 1966b) fall into the first category of dynamic methods. These empirical rela- 
tionsrips were derived on the basis of a crude approximation to the integral over a 
group of plane seismic waves passing by a station. The Gutenberg-Richter estimates 
of energy from agree fairly we- with the static estimates mentioned above. This 
might be exoected, as M s correlates quite well with \og l0 M Q (Kanamori and Anderson, 


1976). 



In this study, we develop dynamic energy estimates of the second kind. We 
apply the theory of Haskell (1964) to compute the energies of several shallow 
events (Section 2), using moments and source-time histories obtained in the last 
decade from sophisticated waveform modeling. Since there are fewer studies avail- 
)!e on intermediate and deep focus events, we also develop a simplified modeling 
procedure (Section 3) to obtain moments and time functions for such events, and use 
these to estimate energy in the same way as for shallow earthquakes. The energy 
estimates we obtain are in a sense direct physical dynamic estimates, as opposed to 
the more empirical approach represented by the energy-magnitude relations, in Sec- 
tion 4, we ct '‘are dynamic and statij estimates for both shaliow and deep events. 

Our dynamic estimates contain more high frequency information than the static 
ones. Thuy are still made, however, at teleseismic distances, and they are further- 
more derived from long period instruments unable to resolve displacement components 
of frequency greater than 1-2 Hz. It is thus possible that some critical high fre- 
quency information is missing. We address this question in Section 6., using high fre- 
quency records obtained close to seismic sources with strong-motion instruments. 

Finally, in Section 6, we compare our dynamic energy estimates with estimates 
from the Gutenberg-Richter energy-magnitude relations, using M s for the shal ow 
earthquakes and long period body wave magnitude m B for the deep and intermedi- 
ate ones. 


2. Dynamic Energy from Source Time Function 

A milestone in the understanding of energy radiation from earthquakes was the 
paper by hcskell (1964). We essentially follow his treatment, with minor modifica- 
tions, to obtain expressions for energy release in terms of parameters obtainable 
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from body wave modeling of earthquakes. The important parameters are the seismic 
moment and the duration and shape of the far field source time function. The earth- 
quake displacement observed at far field is given by 


u(r,f) 


4rrpt/ 3 r 


m o nt) 


(D 


where R(6,tp) is a geometric factor accounting for the radiation pattern of the 
seismic wav'as; p, v, and r are respectively density, elastic wave velocity, and dis- 
tance to the source; Mq is the seismic moment, and T(t ) is the far field source time 
function, which is normalised to unit area. This expression assumes that we have 
already accounted for the effects of attenuation, instrument, receiver structure, and 
geometric spreading (e.g. Langston and Helmberger, 1975). In the simple case of a 
one dimensional rupture with a ramp function near-source dislocation history, T will 
generally be trapezoidal in shape (with a triangle as a special case). The trapezoid is 
obtained by convolving the point-source boxcar (which the near-field ramp produces 
at far field) with another boxcar representing source finiteness. Other shapes are 
certainly possible, though not always resolvable by the data. To calculate the 
energy associated with (1), we begin with a general form of Haskell's (1 964) equa- 
tions (15) and (16) 


• 2ir n 

E = pvj J Ju 2 dt r 2 sind dd dp (2) 

-.0 0 

Equation (2) was derived in the case of spherically symmetric radiation by Yoshiy -na 
(1 963). Rudnicki and Freund (1981) derive it for a more general radiation pattern by 
Imposing plane wave conditions at far field. We apply equation (2) separately to P 
and S waves. We use (1), with R(d,p) factors appropriate (Haskell, 1964) for a 
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OF K- ^ -s'. 

double couple source, and work the geometric integrals out analytically; adding the P 
and S wave energies together, we then obtain 

E = K M§ I t (3) 


WitH 


K = 


1 

1 6npa 3 


+ 


1 

1 0nfi/3 5 


and 


It - ftHt ) df 

where a and /3 are the compressional and shear wave velocities. In the earth, 
/3<a — - so that the 3econd term in K is dominant, and the total energy Is approxi- 
mately equal to shear wave energy. We note that following Plancherel's theorem 
(Bracewell.l 978), (3) can be written as 

E = K Ml I } (4) 


where 


^ = 2/n/) z d/ 

0 

and T (/ ) Is the Fourier transform of T(t) (note that T is real). 

Consider now a simple symmetric trapezoidal far field time function with a ratio 
of rise time to total duration represented hv x (Fig. 1). In this case, the integral ir 


(3) reduct o to 
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/< 


2 

x(i-x) 8 r£ 


( 5 ) 


where T 0 is total duration Hence we have the important result that energy Is pro- 
portional to the square of the moment, and inversely proportional to the cube of the 

duration. If one examines the function — one can easily 'ee that the effect 

X‘1 -x r 

of x Is not Important unless x Is very small; that Is, trape 2 oidal time functions with x 
between ^0.1 and 0.6 have roughly the same energy (Fig 2). When functions have 
very short rise times, this corresponding to the presence of higher frequency com- 
ponents, an appreciable error In the energy can be Incurred from even small errors in 
the rise time. Extremely short rise times are not, however, generally supported by 
the data, and simple but convincing scaling arguments (Kanamori,1 972; Geller,1976) 
lead one to expect values of x greater than 0.1 or so. Hence we effectively have 
two Important parameters In the energy calculation — the total moment end the total 
curation. We might note here that the rather artificial presence of sharp corners In 
the trape/oldal time function does not have an Important effect on the total energy. 
The comers arise from the assumption of a one dimensional rupture. A fault rupturing 
along Its width as well as Its length can be modeled by convolving the point-source 
far-fleld boxcar with two boxcars representing finiteness instead of one, this leading 
to a far field time ‘unction with rounded corners (e g., Mlkumo, 1971, Fig. 2). The 
main shape effect is still due to the rise time, and the above arguments apply. 

We may use (3) to calculate the energies of some shallow events for which time 
functions and moments have been published. Table 1 shows the results of such cal- 


culations, which will be discussed in more detail in Section 4. 
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Figure 


Trapezoidal far field displacement time ' notion. Total duration is Tq, rise 
time is tT 0 . 
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Figure 2. Effect of trapezoid rise time on calculation of dynamic energy release (see 
equation 5). As Song as z (rise time divided by total duration) is greater 
than 0. 1 or so, the effect is not important. 
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TABLE 1 

Energy Calculations for Some Modeled Shallow Events 


Dm* 

<*1*1 

r. 

<wg> 

*• 


1975 

24J 

3 

19.7 


Langston and Butler, 
1976 

1966 

94.8 

3 

19.7 

6.9 

Burdick, 1977 

6/16/76 

25-5 

4J5 

20.6 

6.5 

Cipar, 1960 

9/15/76 
9 21 

24.7 

4.0 

19.2 

6.0 

Cipar, 1980 

9 15/76 
3 15 

25.0 

3-6 

199 

69 

Cipar, 1961 

1967 

25.5 

6.4 

209 

6.4 

Langston, 1976 

1966 

26.7 

4 

219 

69 

Ebel H al , 1978 

1968 

26.0 

6 

219 

69 

Burdick and MeUman, 
1976 

1965 

26 2 

3 

22.7 


Langstcn and Blum, 1977 

1976 

2&2 

6 

21.4 

7.0 

Hartxell, 1980 

1975 

26.6 

7 

22.0 

7.4 

Cipar. 1979 

1975 

27.1 

10 

22.7 

7.7 

Lay and Kanamori, 1980 

7/14/71 

28.1 

14 

249 

79 

Lay and Kanamori, 1980 

7/26/71 

26^ 

16 

24.4 

79 

Lay and Kanamori, 1980 

4/16/65 

25.1 

3.4 

20.5 


Liu and Kanamori, 1980 

9/4/63 

25J2 

Ub 

219 


Liu and Kanamori, 1980 

10/23/64 

25.8 


22.3 


Liu and Kanamori. 1980 

9/30/71 

24.9 

1.6 

21.0 


Liu and Kanamon, 1980 

S/24/70 

252 

Zb 

21.0 


Uu and Kanamori, 1980 

11/29/78 

27.3 

15 

22.6 

79 

Stewart ef al., 1981 

8/23/65 

27.3 

16 

HJb 

7.6 

Chael and Stewart, 1982 

8/2/68 

26.9 

16 

21.7 

7.1 

Chael and Stewart, 1962 

3/14/79 

27.0 

17 

22.7 

7.6 

Chari and Stewart, 1982 

3/24/70 

tbJb 

3 

21.1 

69 

Stewart and Heimberger, 
1981 

1967 

26J 

17 

30.5 

69 

Kanamori and Stewart, 
1976 

1974 

26.7 

22 

20.9 

69 

Kanamori and Stewart, 
1976 


Gibba 
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3. A Simplified Procedure for Modeling Deep Focue Events 

Waveform modeling can be an extremely time-consuming task; the data shown in 
Table 1 represent a very large amount of work on the part of many investigators. To 
obtain a larger data base one may resort to a more simplified procedure which Is still 
sufficiently accurate for the purposes of energy computation. The procedure we use 
Is applicable to deep and intermediate events with comparatively simple sources. It 
consists essentially of estimating the duration of the time function of a simple source 
from the average pulse width of long period WWSSN vertical P waves (Bol- 
linger, 1Q68; Chung and Kanamori,1 980), and then using the average amplitude to 
Infer the moment. We use several stations 10), as well distributed as possible, to 
ave ige out the effects of radiation pattern and directivity. When the long period P 
wave is a single pulse and there are no contaminating free-surface phases, this 
method can be quite accurate. When we applied it to the deep and intermediate 
events studied by Chung and Kanamori (1980), our results for moment and time func- 
tion were in good agreement with theirs. 

To estimate the moment and duration, we use curves of the type shown in Figs. 
3 rnd 4 (see captions). These are obtained from synthetic seismograms which are 
generated by convolving source functions with an instrument response and an 
attenuation filter. We generally assume that the time function is a trapezoid with 
x=0.2 (as we have seen, such a trapezoid does not have a significantly different 
energy from that of a triangle or any trapezoid with x^O.1 ), and T * = 0.7 in the 
attenuation filter. Optimistically, this method, allowing for differences in time function 
shape, attenuation, etc., can give us an estimate of total duration accurate to - 20 
per cent, and an estimate of moment perhaps accurate to within a factor of two, 
given the scatter in amplitude due to receiver and other effects. T he energy 
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Figure 3. (Adapted from Chung and Itanamori, 1980) The relation between measured 
pulse width W p of direct vertical P waves on a long period WWSSN seismo- 
gram and the duration Tq of the far field source displacement time func- 
tion. The curves are obtained by convolving the time function with an 
Instrument response and an appropriate Q filter (7 # * 0.75 shown here). 
These curves are reliable provided the P arrival is a single pulse (ie, the 
event is simple). In this case the event is assumed to be deep enough 
that the direct P wave is not contaminated by free-surface phases. 
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Figure 4. Examples of curves from which the moment M Q can be determined for a 
simple event once the duration r 0 of the far field time function has bae~ 
determined. For a source depth of 400 km., a source-station distance 
60° t and a peak instrument gain factor of 1500, a curve on this diagram 
shows the variation of amplitude Ap of direct P on a long period seismo- 
gram with duration of the time function If the moment of the event 1G* b 
dyn-cm. Thus for a given 7 0 one -an read off the expected amplitude for 
M 0 = 10 es dyn-cm, and compare this with the avetage of amplitude meas- 
urements actually made to obtain the moment of the event (corrections are 
easily made to the amplitude measurement to standardise it to a distance 
of 60° If necessary). Since an average amplitude measurement Is used, 
the curves drawn here are for an average value of the radiation pattern. 
The trapezoid function referred to in the figure has a rise time equal to 
1/6 its total duration, which Is what we generally assume for events we 
are studying by this method. The curves drawn for the limiting cases of a 
boxcar and a triangle show what errors might bt incurred if this assump- 
tion is unwarranted. As can be seen, these errors, as well as those due to 
uncertainties in attenuation, are probably quite negligible compared to 
errors due to scatter In amplitudes caused by receiver and other effects. 
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TABLE 2 

Energies Calculated for Intermediate and Deep Focus Events Studied by Means op 

Simplified Procedure 




Onpi Tb« 



Dspth 

ta) 

■la 

kf M, 
Myna-cm) 

T. 

iaae) 

kf E 
tan) 

M 

D 

Y 

HMjo 

6k 

03 

11 

68 

0626 

32.8 

Tonga-Kermadec 

112 

62 

25.9 

4.6 

20.7- 

06 

12 

67 

0939 

44.3 

Tonga-Kermadec 

134 

6.5 

26.1 

3.7 

21.4* 

12 

06 

65 

1805 

25.2 

Tongs -Kennadec 

156 

6.0 

25.4 

32 

20.0* 

06 

01 

69 

1905 

24.5 

Tongs- Kermadec 

205 

6.1 

25.4 

2.7 

20 2* 

03 

16 

66 

1805 

25.2 

Tonga- Kermadec 

219 

6.0 

25.6 

4.7 

20.1* 

09 

04 

67 

0351 

58.9 

Tonga- Kennadec 

231 

6.2 

25.8 

12 

22.1* 

09 

26 

66 

1437 

46.2 

Tonga • Kennadec 

251 

6.0 

25.3 

2.0 

20.6* 

06 

04 

74 

0414 

13.8 

Tonga- Kennadec 

256 

6.3 

264 

4.7 

21.7* 

02 

22 

75 

2204 

33,5 

Tonga-Kermadec 

333 

6.6 

26.5 

4.6 

21.9* 

01 

20 

68 

2121 

31.6 

Tongs- Kermadec 

349 

6.0 

25.6 

12 

21 2* 

07 

21 

73 

04 'j 

13.7 

T onga- Kennadec 

373 

6.1 

25.8 

22 

21.1* 

06 

27 

70 

i;j5 

18.3 

Bonin la. 

406 

6.6 

27.0 

5.9 

22.4 

11 

16 

66 

2000 

19.5 

Tonga-Kermadec 

424 

6.2 

25.6 

1.75 

21.1* 

11 

29 

74 

2206 

23.5 

Japan 

429 

6.5 

26.6 

6.1 

21.6 

02 

03 

76 

1227 

30.1 

Tonga -Kennadec 

477 

6.0 

252 

2.9 

20.7* 

03 

23 

74 

1428 

33.0 

Tonga -Kennadec 

504 

62 

26.6 

4.95 

21.6 

01 

29 

71 

2156 

03.2 

Japan 

615 

6.6 

26.8 

4.45 

222 

12 

26 

73 

0531 

03.8 

Tonga-Kermadec 

517 

6.5 

262 

2.7 

21 2 

10 

25 

73 

1408 

68.5 

S. America 

617 

62 

25.9 

22 

212 

10 

07 

68 

1920 

20.8 

Japan 

618 

€.7 

272 

13.0 

23.4 

01 

26 

66 

0436 

45.3 

Tonga-Kermadec 

545 

5.8 

252 

L75 

202* 

01 

24 

69 

0233 

03.4 

Tonga-Kermadec 

587 

6.7 

26.1 

0.75 

23! 

06 

26 

70 

1109 

51.3 

Tongs -Kennadec 

587 

6.1 

25.7 

2.35 

20.9 

07 

21 

66 

1830 

15.3 

Tongs -Kennadec 

590 

5.8 

25.8 

12 

212* 

02 

16 

67 

1011 

11.8 

S. America 

598 

6.4 

262 

4.1 

212 

10 

09 

67 

1721 


Tonga-Kermadec 

805 

6.8 

27.0 

4.9 

22.4 

03 

24 

67 

0900 

20.0 

Java 

601 

62 

26.1 

4.1 

20.9 

03 

17 

66 

1650 

33.1 

Tonga-Kermadec 

630 

62 

262 

4.0 

21.6* 

10 

01 

72 

2349 


Philippines 

632 

52 

25.0 

0.7 

21.0 

02 

10 

60 

2258 

03.3 

Tonga-Kermadec 

635 

6.4 

25.6 

52 

21.7 

12 

09 

65 

1312 

65.3 

Tonga-Kermadec 

649 

5.7 

25.7 

12 

212* 


* Events studied by Chung and Kanamori (1960). 



- 19 - 


0R1G1NAL PAGE 
OF POOR Q^ A - 


r ,G 

iTY 


estimate is probably good to an order of magnitude or so. Energies calculated for 
deep and intermediate events studied by this method, including the events of Chung 
and Kanamori (1 980), are listed in Table 2. 


4. Comparison with Static Energy Estimates 

We now examine the results of the energy calculations in the framework of an 
important independent method of estimating energy, based on elastostatic considera- 
tions. Consider a simple model of an earthquake where a 0 , cr lt and are initial, final, 
and dynamic frictional stresses on the fault. We may write (Savage and Wood, 
1971) 


W = 


CT 0 + g i 



D S 


( 6 ) 


where W is the difference between the strain energy drop and the frictional energy, 
D is the average dislocation, and 5 is the slip area. By using the stress drop 
Act = n 0 -ci and the seismic moment A/ 0 = i±DS, we can rewrite (6) as 


W = 


A a t ( g i~ g /) 

2/i /i 


M o 


(7) 


Orowan (1960) proposed a physically very reasonable model of a fault whereby 
motion stops when the accelerating stress decreases to a value equal to some aver- 
age dynamic frictional stress, i.e. a, = t ij. There is thus no overshoot arising from, 
say, the inertia of the moving fault blocks. In Orowan's model eq. (6), which is the 
strain energy drop less the frictional energy represents the energy radiated as 
seismic waves. If Orowan's condition is satisfied, then clearly the second term in (7) 


vanishes, and we have simply 



ORIGtJV.L . v; 
OF POOR v- ..." 


w = 


A crM 0 

~r;r 


( 8 ) 


Kanamori (1977) used this relationship to estimate the energy released in great shal- 
low earthquakes. With Aa » 20 to 60 bars (2 to 6 X 10 7 dyne/cm 2 ), and 3 to 6 
X 1 0 1 1 dyne/cm 2 ), 


r 0 «(5A'1O- 8 >A/ 0 (9) 

where we have now adopted the subscript 0 to indicate that this is a static or 
essentially ze r o frequency estimate of energy, as opposed to the higher frequency 
estimates made from (3). 

Fig. 5 shows a plot of energy vs. moment for the shallow events of Table 1. The 
line shows the energy according to (9), with parallel lines bounding an order of magni- 
tude up or down. There is considerable scatter. Some of this scatter must be due to 
the errors in To and M 0 ■ Another contributing factor, however, probably arises from 
the fact that (9) is derived assuming Act is 20-60 bars, and for many events this 
obviously need not be true. The dynamic estimates by their very nature take into 
account the details of rupture for the individual events. For this reason, they can 
deviate considerably from the line E = (6X1 0~ 5 *Af c , perhaps even more than would a 
crude estimate from M m . An Interesting case is that of the two Gibbs fracture zone 
events (Kanamori and Stewart, 1976). They lie considerably below the line. As they 
are known to have been especially slow events, it should not be surprising that (9) 
might overestimate their energy. 

All in all, considering the simplicity of the model leading to the static estimate, 
the errors in the dynamic estimate arising from errors in M 0 and T c , an.' the indepen- 
dence of the two methods, the agreement between the static and dynamic energy 
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Figure 6. Energy calculated for some modelled shallow events (Table 1) plotted 
against seismic moment. The line shown corresponds to the approximate 
relation F - (6X10 S \W 0 (which assumes a stress drop of 20-60 bars) 
obtained by Kanamori (19 77). The parallel lines bound an order of magni- 
tude up or down. Considering that this simple elastostatlc calculation Is 
completely Independent of the dynamic calculations made here from body 
waves, the agreement Is encouraging (see Section 4). 
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determinations for shallow earthquakes is rather good. We may examine this rough 
equality more closely by considering some simple static stress drop scaling relations. 
In the case of constant stress drop, we may write the moment in terms of stress drop 
and fault area as (Kanamori and Anderson, 1975) 

«0 = A aS 3/z (10) 

Using an approximate expression ro w “ — for the time function duration, we obtain 

P 

( 11 ) 


Substituting this into (3), and using (5), we have 


E * 


2K 

x(1 -z) 2 


A o/PMo 


( 12 ) 


Using x = 0.2, /? = 3.4 km/sec, A a = 30 bars, and p * 2.8 g/cc in K gives us 

E a (4.6*1 O- 5 )M 0 (13) 


which is very close to (9). 

Fig. 6 shows energy versus noment for the deep and intermediate events listed 
in Table 2. The lines are the sat.ie as the ones in Fig. 5. On the whole, the deep 
events tend to plot below the line corresponding to W D = (5X'\0 ~ S )M 0 Of course, 
given that our energies are not likely to be accurate to better than an order of mag- 
nitude, this may not be significant. However, the effect is quite systematic, and con- 
trary to what one might expect if one believed that deep events tend to have higher 
stress drops: the average stress drop determined by Chung and Kanamori (1980/ for 
their deep and intermediate events is ^ 500 bar. If p. * 6 to 10 X 10" dyne/cm z 
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Figure 6. Similar to fig. 6, but for deep and intermediate events. Circles are for 
events In Table 2; open circles in particular are for events also studied by 
Chung and Kanamori (1980), and closed ones are for the rest. Closed 
squares = Mikumo (1971). Closed triangle * Fukao (1970). 
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below 400 km, the relaton W ' * (6 X 1 0 s ) A/ c would require Aa ^ 60 to 1 00 bars, so 
If one believed the high stress drops of Chung and Kanamori (1980), one would 
expect at least the events they studied (we have not determined stress drop for the 
extra events we studied) to plot above the line. 

The key to understanding this situation may lie In remembering that for (8) to 
hold, Orowan’s condition must be met, and this need not be the case. If we assume 
that the condition is met, we may solve (8) for Aa, and use values of moment and 
dynamic energy to obtain a value of stress drop which we may call 'Orowan stress 
drop*. This value should be equal to the actual stress drop If Orowan’s condition Is 
met; If not, It should be lower. If we calculate Orowan stress drops for the events of 
Chung and Kanamorl (1980), we find that they are considerably lower (Fig. 7) than 
Chung and Kanamorl’s teleseismically calculated stress drops (using inferences of 
fault area from the source time functions). If we calculate the Orowan stress drops 
using energy determined from m# (see Section 8) Instead of our dynamic estimates 
from Section 2, the gap is even wider. The implication, then, Is that either Orowan's 
condition is not met for these events, or the condition Is met and the Chung-Kanamori 
stress drops are too high, by almost an order of magnitude. Since stress drop Is one 
of the more model-dependent and poorly determined seismological quantities, this 
would not be too surprising. 

In any case, It Is not difficult to see why a relationship of the form 

t: = qU c (14) 

can hold for deep and shallow events alike with q approximately given by 5 X 10 ! \ 


From ( ?> we see that 
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Flgure 7. Upper dotted line * Stress drops obtained by Chung and Kanamorl (1980) 
(using source dimensions Inferred from far field time functions), plotted 
against depth. Lower dotted line * 'Orowan stress drops', calculated from 
equation (8) assuming Orowar.'s condition Is met, and using energies 
obtained from m B . Solid line * 'Orowan stress drops' calculated from 
equation (8) assuming Orowan's condition is met and using dynamic ener- 
gies calculated in this study. The use of our energies , which are generally 
higher than those estimated from m B , does not rlose the gap between the 
Orowan stress drops calculated from (8) and those obtained by Chung and 
Kanamori (1980). Either (1) our energies are systematicaly too low or (2) 
Chung and Kanamori's (1980) stress drops are too high or (3) Orowan's 
condition is not met for these deep and intermediate events. 
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9 = 



M 


(15) 


in the case of shallow events, where Orowan's condition is likely to be met (Kanamori 
and Anderson, 1975), we merely have the reasonable condition, as stated before, 

that 6X1 0~°. For deep events, we can have a similar situation as for shallow 

events, or we can have a non- Orowan process with high stress drops In the first term 
of (15), and a negative second term. 


6. Near Source Energy Studies and the Question of Frequency Content 

The computations which we have carried out are based on earthquake displace- 
ment data viewed through a variety c f distorting filters, such as attenuation and 
instrument. We address here the question of the validity of these results, given that 
by using a long period instrument cannot hope to resolve displacement com- 
ponents of frequency greater than 1 to 2 Hz. Beyond the problem of the instrument, 
we must also consider the possibility that important high frequency energy is 
attenuated, either anelastically or through scattering, by propagation to teleseismic 
distances. One could make the argument that high frequencies observable only very 
close to the source could be responsible for a considerable portion of the total 
energy. We note here that we cannot simply quote the fact that teleseismic comer 
frequencies are relatively low for earthquakes of size similar to the ones examined 
here as evidence that high frequencies are unimportant. A teleseismic spectrum is 
not necessarily simply related to the true source spectrum at near field. 

An important source of information with regard to these questions is to be found 
In near-source strong-motion records. By examining data obtained c'ose (< 20 km.) 
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\o the source using high frequency strong-motion Instruments, we can assess the 
Importance of the shorter period energy. From an accelerogram, one can easily obtain 
a velocity trac and use that to compute the quantity 


/)(/) = / u' df ' 


( 10 ) 


which is proportional to the integral of the energy spectrum to a given frequency; 
u(f) Is tho Fourier transform of the velocity trace The seismic wave energy 
obtained from a trace at a given station is given approximately by 

E * Anpftr* R (6,if) 2 D(°°) (17) 


We note that (16) is not the integral of the source energy density per se, but of 
the trace energy density. We are thus not looking directly at ts.o true source spec- 
trum. There is some contamination from reflection, refraction, scattering, etc. How- 
ever, if the high frequency contribution in traced obtained close to the source is not 
Important, i.e. if D at 2 hz appears to have already reached a final value, then we can 
probably not bo too worried that we are looking at a trace spectrum rather than a 
tr je source spectrum. That is to say, If large amounts of high froquency energy were 
present, we might have to be concerned that the contaminating processes we have 
mentioned might be the origin of It, but if such energy is not there It does not matter 
as much to our argumenv that such processes might be present. The contaminating 
processes we have mentioned would probably, If anything, enhance the high fre- 
quency content of the trace relative *o the source, which by Itselr would argue that 
If high frequency energy is negligible in the trace, It must also be negligible in the 


sou ce. Of course, this ignores attenuation; if we are close enough to the source, 
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however, attenuation should not be Important. We discuss this more fully below. 

Figs. 8 a,b,c,d show [j(f ) for several records from the 1 97 1 San Fernando, 1 966 
Parkfield, 1979 Imperial Valley, and 1977 Bucharest earthquakes. Table 3 shows 

values of D( 10) and the ratios » J j(fo y~ an< * p^Q) ' w ^ 8re * he argument is in 

Hz, for these and other records. We use D( 10) to be essentially representative of 
D( no). This certainly seems justified on Inspection of the figures (In addition, sam- 
pling Intervals for the digital data are often such that folding frequencies themselves 
are not much higher than 10 Hz. ). Many of the records were obtained extremely 
close to the source (e.g. Pacolma, less than 1km. from the nearest point on the 
Sierra Madre Fault (Heaton, 1932)), and in r o case is any apr reciable energy 
observable above 4 Hz. Such energy may exist in the very immediate vicinity of the 
source, but In that case we may raise semantic questions about which energy to con- 
sider "radiated" and which not. If this hypothetical high frequency energy is 
attenuated within 1 km. of the source, we cannot consider it to be radiated energy. 
This reasoning applies also to energy at 1 to 2 Hz. If there Is Important energy in this 
band which we cannot see even at 1 km or so from the fault (actually, with a Q of 
about 300 this Is unlikely), then we can hardly worry about it for the purposes of 
computing radiated seismic energy. 

What we have set out to do in examining the strong-motion records is to see if 
there was a large proportion of energy there which we were missing at teleseismic 
distances. It is clear from the records presented here that even very close to the 
source, by far the largest proportion of the energy is contained in frequencies below 
2 Hz. In many cases, over 90 per cent of the energy Is even below 1 Hz. What these 
results suggest Is that no appreciable error (certainly not one of an order of magni- 


tude) Is incurred by making an energy determination at far field using a long period 
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Flgur© 8. The Integral D(f) (see section 6) of the spectral energy density versus 
frequency from strong motion velocity (cm/sec) traces tor the San Fer- 
nando, Parkfield, Imperial Valley, and Bucharest earthquakes. Different 
curves for each earthquake correspond to different records, (see Table 3) 
The curves for the San Fernando, Parkfield, and Imperial Valley earth- 
quakes suggest that, even close to the source, by far most of the energy 
radiated Is below 1 to 2 Hz, In frequency. Far field energy determinations 
using long period Instruments thus may not be in gross error, despite their 
bandtimlted nature. 
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instrument. 

Strictly, this only applies to shallow events. Certainly we have no instances of 
strong motion recordings within 1 km. of the source of a deep focus event, so we 
cannot directly address the problem of whether there is important energy within a 
few kilometers of the source which never propagates out to teleseismic distances. 
We can, however, make some statement about whether or not a long period instru- 
ment is bro .d enough in its frequency response to retrieve adequately the energy 
that does manage to propagate tc the teleseismic range. The curves of Fig. 8(d) for 
the 100 km. depth Bucharest earthquake in fact show very little energy outside the 
passband of a long period WWSSN instrument (fc 60 sec. to 1 to 2 Hz.), and this is 
encouraging. 


6. Energy and Magnitude 

In this section we compare our dynamic energy estimates with the energies one 
would obtain using the Gutenberg-Richter relations. For the shallow events of Table 
1, the comparison is relatively straightforward; we may use M s as a measure of mag- 
nitude. Fig. 9(a) shows log 10 £ in ergs versus M s for these events. Our estimates 
seem to be consistently lower than the Gutenberg- Richter line. A best fit line 
through our points would have slope 1.81 (i 0.2) and intercept 9.06(± 1.38), com- 
pared to 1.5 and 1 1.8 respectively for Gutenberg-Richter. 

The comparisor for the deep and intermediate events of Table 2 is more ambigu- 
ous. These events generally did not excite appreciable surface waves, so we must 
use a body wave magnitude. Gutenberg and Richter (1956a,b) derived the relation 
log }0 E = 2.4m£ + 5.8 The magnitude ttlq is not the same as the m b now in common 
use. The latter is a short period (% 1 sec) body wave magnitude, while the former is 
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a longer period one. We have used long period WWSSN records to determine an m g 
more compatible than jn b with Gutenberg and Richter's definition. 

One difficulty which arises is that when the P wave consists essentially of a 
single pulse, as is the general case with the simple events we have studied here, 
the measurement of the dominant period in the wave group becomes ambiguous. We 
have set the perioa io twice the pulse width. Another difficulty is that the WWSSN 
instruments whose records we have employed are peaked at 1 5 seconds, while 
Gutenberg and Richter used mechanical instruments with a different period response 
(flat rather than decaying); thus, one must be careful to use the correct gain for the 
WWSSN instrument when one is looking at a period different from the peak period. 
The waveforms from the two instruments differ; we have conducted some numerical 
experiments to ascertain that no drastic errors occur because of this. 

A plot of log E versus m B for the intermediate and deep events of table 2 is 
shown in Fig 9(b). In contrast to the case of the shallow events, the bias here is 
above the Gutenberg-Richter line. The least squares line through our plotted points 
has slope 1.97(± 0.34) and intercept 9.07(± 2.13). We note that if one allows an 
error of 0.5 units in m B , taking into account all the factors mentioned above, as well 
as an error of an order of magnitude in the energy, the discrepancy is understand- 
able. 

Although it is interesting that the shallow events generally plot below the 
log l0 £’— A/j line, while the deep and intermediate ones plot above the log 10 E-mg 
line, we cannot really make meaningful comments about this given the empirical nature 
of the Gutenberg-Richter relationships. 



40- 


Figure 9. (a) Common logarithm of the dynamic energy release in ergs plotted against 
U, for shallow events of Table 1. The line represents the Gutenberg- 
Richter relationship. 
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Figure 9. (b) Common logarithm of the dynamic energy release In ergs plotted against 
m B (long period body wave magnitude— see section 6) for the deep and 
intermediate events of Table 2. Squares represent events also studied by 
Chung and Kanamori (1980), The iine represents the Gutenberg-Richter 
relationship. 
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7. Conclusion!* 

(1) The important parameters in the calculation of seismic emorgy release from 

body waves are seismic moment fj 0 and far field displacement time function duration 
M 2 

7'n. with E a — The Important shape effect tor the usual trapezoidal time function 

TS 

comes from its ratio z of rise time to total duration. A v . as 0.1, which is gen- 
erally supported by the data, the effect is not important. 

(2) Our near source studies suggest that most of the important radiated energy 
is below 1 to 2 Hz. in frequency, and hence that far field energy determinations using 
long period WWSSN instruments are not in gross error despite their bandlimited 
nature. 

(3) Dynamic energy estimates for shallow earthquakes made from body waves 
are in reasonable agreement with expectations from simple static elastic relaxation 
models, which suggest that E^5X 1O” 5 A/ 0 for shallow events when a stress drop of 
20-60 bars is assumed. 

(4) Deep events, despite their possibly different seismoioglcal character, yield 
dynamic energies which are also '.T'mpatible with a static energy pre ictlon similar to 
that for shallow events. Seismic moment J and hence a moment based magnitude 
scale, may reliably be used for shallow and deep events alike, as a reasonably accu- 
rate measure of energy release. 
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PART II: SUBDUCTION ZONE SEISMICITY AND STRESS IN 

SLABS 
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Abstract 


The pattern of seismicity as a function of depth in the world, and the orientation 
of stress axes of deep and intermediate earthquakes, are explained using viscous 
fluid models of subducting slabs, with a barrier in the mantle at 670 km. 670 km Is 
the depth of a seismic discontinuity, and also the depth below which earthquakes do 
not occur. The barrier in the models can be a viscosity increase of an order of magni- 
tude or more, or a chemical discontinuity where vertical velocity is zero. Log N 
versus depth, where N is the number of earthquakes, shows (1) a linear decrease to 
about 260-300 km depth, (2) a minimum near that oepth, and (3) an increase 
thereafter. Stress magnitude in a subducting slab versus depth, for a wide variety of 
models, shows the same pattern. Since there is some experimental evidence that N 
Is proportional to e kc , where A: Is a constant and a is the stress magnitude, the 
agreement Is encouraging, in addition, the models predict down-dip compression in 
the slab at depths below 400 km. This has been observed in earlier studies of earth- 
quake stress axes, and ve have confirmed it via a survey of events occurring since 
1977 which have been analysed by moment tensor Inversion. At intermediate depths, 
the models predict an approximate but not precise state of down-dip tension when 
the slab is dipping. The observations do not show an unambiguous state of down-dip 
tension at intermediate depths, but in the majority of regions the state of stress Is 
decidedly closer to down-dip tension than it is to down-dip compression. Chemical 
discontinuities above 670 km, or phase transitions with an elevation of the boundary 


In the slab, predict, when incorporated into the models, stress peaks which are not 
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mlrrored In the profile of seismicity versus depth. Models with an asthenosphere and 
mesosphere of appropriate viscosity can not only explain the state of stress 
observed In double Benloff zones, but also yield stress magnitude profiles consistent 
with observed seismicity. Models where a nonlinear rheology Is used are qualitatively 
consistent with the linear models. 


1. Introduction 


In this study we use simple models of subducting slabs to explain observations 
of the distribution of earthquakes versus depth, and observations of the orientation 
of stress axes of deep (> 300 km) and intermediate (70 to 300 km) earthquakes. 

The distribution of earthquakes with depth has been discussed by many investi- 
gators (e.g., Gutenberg and Richter, 1964; Sykes, 1966; Isacks et al, 1968), who 
have variously noted the presence of seismicity minima near 260 km, and the 
existence of deep peaks in seismicity. Richter (1979) has recently explored the 
possibility of a barrier to mantle flow at the 670 km seismic discontinuity (e.g., 
Whitcomb and Anderson, 1968) being responsible for the large increase in seismicity 
above this depth in the Tonga-Kermadec region. He also argues that such a barrier 
may explain the tendency toward down-dip compression at depth for earthquakes in 
this region (previously noted by Isacks and Molnar, 1969, 1971). 

We proceed here in the same spirit as Richter extending the observational base 
to the whole world, and exploring a wide range of models of subducting slabs. In 
section 2 we present observations of seismicity versus depth, using the large amount 
of catalog data that has become available since the studies of Sykes (1956) and 
Isacks et al (1 968). Wn also analyse the orientations of stress axes of earthquakes 
studied using the new methods of moment tensor inversion (Kanamori and Given, 
1981 ; Dziewonski et al, 1981; Appendix A). In Section 3 we present calculations of 


stress in subducting slabs. 
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2. Observations 

2.1 Seismicity as a Function of Depth 

Fig. 1 shows histograms of the logarithm of the total number of earthquakes in 
the world versus depth. The different curves correspond to different cutoff magni- 
tudes: one represents the earthquakes with a o r «e-seeond body wave magnitude m b 
greater than or equal to 4; In the other two, the cutoffs are m 6 = 5 and m b = 6. The 
time period covered is 1964 to 1980, 1964 being the year in which the one-second 
body wave magnitude began to be applied uniformly. The data sources are the NOAA 
(1964-1977 inclusive) and PDE (1978-1980) catalogs. 

The curves have some striking features. First, we note the well established 
fact (Isacks et al., 1968) that there are no earthquakes below a depth of about 700 
km. Second, we observe a roughly linear decrease in log N from the surface to a 
depth of approximately 260 - 300 This exponential behavior was noticed by 
Sykes (1966), and Isacks et al. (1968); it has not been discussed very much in the 
literature since. After the linear- og decrease, there appears to be a seismicity 
minimum, followed by a resurgence in activity from 600 to 700 km. The three curves 
behave very similarly in these respects, although the curve for m b > 6 is spottier 
than the other two, probably showing the effects of incomplete sampling of larger 
earthquakes in the time period covered. 

The worldwide curves of fig. 1 represent what might be termed an ' average 
subduction zone . Although these curves contain very important information, it is 
necessary also to examine similar plots for individual subducting regions. Such plots 
are presented in Appendix C. In deciding how to divide up the regions, we have relied 
on the physiography of trenches, the planform of seismicity observable on global 
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Flgure 1. Logarithm of the total number of earthquakes in the world versus depth, in 
20 km intervals. The sources are the NOAA and PDE catalogs, from 1 964- 
1 980. The three curves are for throe different cutoff magnitudes, as 
noted, where denotes the one-second body wave magnitude reported 
by the ISC. 
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maps such as that of Tarr (1974), the stereo plots of Johnson and Richter (1979), 
and previous studies of the lateral segmentation of subduction zones (e.g., Isacks 
and Barazangi, 1977). We have tried to isolate portions of regions where the struc- 
ture of the Wadatl-Benioff zone is as close to two-dimensional as possible. 

If we examine the seismicity curves of individual regions, we find that many of 
them faithfully reproduce the global pattern in whole or in part. Regions without 
deep seismicity tend to roproduce the upper part of the pattern, that Is, the roughly 
linear-log decrease down to 260-300 km. Regions with deep seismicity tend to 
reproduce the entire pattern, with the roughly linear decrease down to 260-300 km., 
followed by a minimum and a deep peak of varying position and intensity. 

As we ^an see from the plots, however, several regions do not conform precisely 
to these specifications. Many of the shallow regions have a decay pattern of seismi- 
city with depth which deviates significantly from linearity. Note the pattern in the 
Hindu-Kush, which is more a zone of continental convergence than a subduction zone 
(Molnar and Tapponler, 1975). Some of the shallow dipping South American zones 
show patterns more similar to this than to the patterns of other circum-Paciflc zones. 

To first order, however, the features of the seismicity curves which we have 
listed do appear to be global, and bear explanation. 
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2.2 Orientations of Stress Axes of Intermediate and Deep Focus Earthquakes 
2.2.1 A First Look 

We now turn our attention to observational evidence concerning the orientation 
of stress in subducting slabs. In a pair of classic papers, Isacks and Moinar 
(1969,1971) examined the relationship of earthquake stress axes, lerived from first 
motion studies, to the geometry of Benloff zones. In the 1969 paper they established 
the result that stress axes (tension or compression) are more closely aligned with 
slab geometry than are nodal planes. That is, the evidence does not suggest that 
earthquakes represent shear motion along a fault plane marking the Interface 
between slab and mantle. Actually, some evidence for this has been discovered (e.g , 
Um'no and Hasegawa, 1982), but this is for depths shallower than 60 km. 

In the 1971 paper, Isacks and Moinar performed a detailed regional analysis. 
Globally, the evidence seems to suggest down-dip compression or tension at depths 
from 70 to 300 km, depending on the region, and down-dip compression below 300 
km. Of course, ’down-dip ' is to some extent in the eye of the beholder. In this sec- 
tion and the next we will attempt to quantify how close to down-dip the stress axes 
of earthquakes are, in the mean. Fujita and F.ananori (1981), who have performed a 
global survey of focal mechanism soluJons for intermediate earthquakes, and had 
more data at their disposal than did Isacks and Moinar, argue for a state of 'in-p.ate 
rather than down-dip stress. 

First-motion solutions vary greatly in the degree to which the nodal planes and 
stress axes can be constrained by the data. Heaton (1982) has compiled a very 
useful table of the best constrained solutions for earthquakes deeper than 30 km. 


Several of Isacks and Molnar's solutions are in* jded, as are first -mot ion studies by 
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more recent investigators, and some more detailed studies involving the use of syn- 
thetic seismograms, in fig. 2, we plot the stress axes from Heaton's list on lower- 
hemisphere stereographic projections. The data from different repions are combined 
by rotsting all axes such that the slab is always vertical and striking North. The dips 
and strikes of the various slabs are taken from Uyeda and Kanamori (1979) and 
Fujita a.id Kanamori (1981). There is some scatter introduced by the fact that slabs 
are not everywhere even approximately two dimensional 'n structu-T, but we have 
tried to minimise this by eliminating events where the average dips and strikes can- 
not clearly be used. Particular care must be taken in the Tonga, Indonesian, and Phi- 
lippine regions. 

Looking at fig. 2, we see considerable scatter In the data. However, a general 
tendency towards down-dip compression is definitely observable in the deep events. 
In the Intermediate events there is less of an obvious cluster, igh the tension 
axes do S'’em to line up with the down-dip direction better than does anything else. 
If we choose for erch earthquake the stress axis closest to the down-dip direction 
and plot that axis, then the cluster improves for intermediate events, altho'-gh not for 
deep ones. This is because of the more universal nature of down-dip compression at 
depth, as opposed to the mere region-dependent state of stress above 300 km. In 
any case, we can see that the stress axes, rather than the nodal planes, are what 
align best with the geometry of the slab. 
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Figure 2. Stereographic (equal-area, lower hemisphere) projections of stress axes 
and nodal plane poles for earthquakes drawn from Heaton's (1982) list of 
well constrained solutions. Some of (sacks and Molnar's (1971) mechan- 
isms are included, as well as more recent onas. Figure 2(a) is for deep 
300 km) earthquakes, and figure 2(b) is for Intermediate (70 to 300 km) 
earthquakes. Ail quantities are plotted in a slab coordinate system. The 
slab Is always vertical and striking North-South. The down-dip direction is 
at the center of the stereonet. Compression and tension axes are shown. 
Also shown is a drawing where the stress axis (compression or tension) 
closest to the down-dip direction is plotted. 
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2.2.2 Events Studied by Moment Tensor Inversion 

An important new development in seismology has been the development of 
methods for performing routine moment tensor inversions for earthquakes using now 
available digital data. In the years preceding the availability of this data and the 
development of the new methods, an Invostigator who wanted to do a more detailed 
study of an earthquake mechanism than a first-motion analysis had to commit a size- 
able amount of time to the digitization and forward modelling of the records. The 
number of events that couiu be studied in this manner was therefore quite limited. 
The new moment tensor inversion methods offer the advantages of objectivity, as 
well as thoroughness and speed. We have studied several of the larger (A/q^IO 20 
dyr.-cm) intermediate and deep events which have occurred in the world since 1978, 
using the method of Kanamori and Given (1981) for inverting IDA data. Details are 
given in Appendix A. In addition, we have available several solutions performed in 
independent studies by Dziewonski and Woodhouse (1983) and Giardini 
(1962,1983). These investigators have obtained reasonable solutions for events 
with M o&1 0 24 dyn-cm. Their lower threshold stems fror their effective use of body 
w&ves from the SRO; our IDA inversions are carried out on surface waves at periods 
exceeding 200 seconds. We now ask ourselves: do these earthquakes studied in 
this more objective and thorough manner lend support to the general conclusions dis- 
cussed in the last section? 

In general, they appear to. Fig. 3 shows the stress axes for intermediate and 
deep events studied by moment tensor inversion, plotted as before, with data from 
different regions rotated such that the slab is striking North and vertical. In this set 
of figures, too, the most readily apparent feature in the considerably scattered data 
Is the tendency toward down-dip compression for the deep events. In the 
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Figure 3. As Figure 2, but for events from 1977-1681 analysed by moment tensor 
Inversion, in this study (Appendix A), and in Dziewonski and Woodhouse 
(1983) and Giardini (1983). Figure 3(a) is for deep events, and Figure 
3(b) is for intermediate events. We recall that down-dip is always in the 
center of the stereonet. 
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Intermediate events, a pattern is more difficult to discern, but the tension axes do 
show the tightest group around the plate in general, If net convincingly around the 
down-dip direction. 

Fig. 4 shows the results of some simple statistical analysis of the data in fig. 3. 
Since we are dealing here with the statistics of directions rather than conventional 
linear statistics, some special techniques must be applied. The techniques are well 
developed and have long been in use in paleomagnetism and biometry; a brief 
description and references are given in Appendix B. Fig. 4 shows both Bingham and 
Fisher statistics for the tension and compression axes of tha intermediate and deep 
focus earthquakes. The Bingham and Fisher means are quite similar. In all cases, the 
hypothesis of uniformity can be rejected to 99% level or better. This means that 
there is 1% chance or less that the data are drawn from an isotropic distribution, and 
a preferred direction does not exist. The larger circle or ellipse in each piot 
represents the standard deviation of the data, while the smaller one represents the 
boundary of the region of 95% confidence for the mean direction. 

From these figures we can conclude, despite large s itter, that the deep 
events are consistent with down-dip compression but not down-dip tension, while the 
intermadiate events are marginally consistent with down-dip or in-plate tension, but 
not down-dip compression. 

We must remember, however, that our sample contains data from many different 
regions. While we do not have enough events to conduct an exhaustive individual 
survey of each of the world's subduction zones, it is still worthwhile to present some 
of the data by region, or class of region. One important subdivision we must perform 
for the deep earthquakes is to isolate the Tonga region. This is the most active area 
of the world at depths below 300 km and accounts for 46% of the deep events in 
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Flgura 4. Some simple statistics! parameters for the data shown in figure 3. Both 
Fisher and Bingham statistics are shown (Appendix B). For convenience of 
plotting, th*> projections are now Wulff projections, and cover the whole 
sphere. An open square Indicates a mean position in the upper hemisphere, 
while a filed square Indicates one in the lower hemisphere. Small ellipses 
or circles show 95% confidence limits for the true mean directions. In all 
cas* a preferred or mean direction does exist, to 99% confidence or 
better. Large ellipses or circles delineate standard deviation limits. 
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our sample. Fig. 5 plots stress axes for deep earthquakes separately for Tonga and 
for the rest of the world. The tendency toward down-dip compression at depth is not 
restricted to Tonga, as we can see from the statistics In fig. 6. Fig. 7 shows 
compression axes for individual regions where we have more than four events. The 
number of events in each case is small enough such that it is difficult to draw a con- 
fident conclusion for the region. We can see, however, that North Honshu, Mindanao, 
and Izu-Bonin are all more or less consistent with down-dip compression. The large 
scatter in the case of Izu-Bonin is probably at least partially due to the assignment 
of an average dip to the entire subduction zcne, whereas this zone appears to 
change dip from North tn South (Katsumata and Sykes, 1969). This may also be a 
factor in the case of Java, whose state of stress is not well resolved in the figure. 

Turning now to intermediate depths, we are faced with the problem that Tonga 
and the New Hebrides are really the only regions which ar r individually well 
represented in our sample. Tonga shows some tendency toward down-dip compres- 
sion, in agreement with the findings of Richter (1979) and Isacks and Molnar (1971) 
(see fig. 8). The New Hebrides show some tendency toward down-dip tension, as 
was also observed by Pascal et al (1978), as well as Isacks and Molnar. In both 
cases, actually, "in-plate” might be a more accurate expression than "down-dip ' 
(Fujita and Kanamon, ISB 1 ). Shown also are the axes for North Honshu, the Maria- 
nas, and Mindanao. These are all deep extending slabs, but they do not show down- 
dip compression at intermediate deoths the way Tonga dce^. Mindanao seems closer 
to tension than compression (see also Cardwell et al. , 1 980). However, these 
regions are not well represented at a!i in our sample, and we cannot make any strong 
statement about them. Other regions are represented, individually, even worse. 
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Figure 5. As figure 2, but plotting axes seoarately for Tonga and the rest of the 
world. 
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Figure 6. Bingham statistics, plotted as fn figure 4, for the data representing the 
world except for Tonga (figure 5). 
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One potentially instructive way to look at the data is shown in figs. 9 and 1 1. 
Here we plot stress axes for intermediate depth earthquakes occurring in slabs that 
do not extend below 300 km (fig. 9), and for intermediate depth earthquakes occur- 
ring in deep-extending slabs other than Tonga (fig. 10). Again, the scatter is large, 
and it is difficult for us to draw conclusions as confidently as we have for deep 
earthquakes As we can see from the mean directions in fig. 10, shallow-extending 
slabs are closer to down-dip tension than compression. This also appears to be true 
of the deep-extending slabs other than Tonga (fig. 12). This disagrees with one of 
the conclusions of Isacks and Molnar (1971), who believed that deep extending 


slabs in general were in down-dip compression at all depths. 
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Figure 7. As figure 2, for deep events in individual regions where more than 4 events 
are available. 
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Flgure 8. As figure 2, for Intermediate events in Individual regions where more than 4 
•vents are available. 
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Figure 9. As figure 2, for Intermediate events in slabs with maximum depths at 300 
km or less. 
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Flgura 1 0. Bingham statistics, plottad as in flgura 4, for avants In shallow-axtandlng 
slabs (flgura 9). 




Figure 11. As figure 2, for Intermediate earthquakes In slabs other then Tonga with 
maximum depths below 300 km. We separate Tonga because of its ten- 
dency, visible in figure 8, toward down-dip compression. 
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Figure 1 2. Bingham statistics, plotted as In figure 4, for data shown in figure 1 1. 





3. Calculations of Stress In Subducting Slabs 


3.1 Introduction to tho Models 

We presume that, other things being equal, the number of earthquakes which will 
occur In a given area is related to the level of stress. Section 3.4 discusses this In 
more detail. We thus adopt some simple models of a subducting slab and calculate the 
stresses therein. We must emphasize that what we seek is not a detailed numerical 
simulation of any individual subduction zone — a futile tack given the complexity and 
many unknowns involved— but a series of models which can elucidate the basic phy- 
sical processes and give us qualitative Insight. 

It has long been thought that at long time scales the Earth's mantle behaves as 
a fluid, and that plate tectonic processes are associated with a large scale thermally 
driven mantle circulation. The literature on this topic is extensive. Hager and 
O'Connell (1981), O'Connell (1977), and McKenzie et al. (1974) are good starting 
points for one who is Interested. We will not be reviewing here the solid state 
processes by which a solid nantle might creep or flow over geologic time. The litera- 
ture here Is also extensive, but we can quote Ashby and Verrall (1977) and Gueguen 
and Nicolas (1980) as general references. 

Our models consist of a box of fluid, as shown in fig. 1 3. The slab is modelled as 
a denier and more viscous fluid than the surrounding mantle, as we will discuss 
shortly. We solve the Stokes problem for viscous Incompressible flow via a penalty 
function, finite element method (lii'ghes et al., 1979). Typically we use grids of 
square or rectangular elements. We always have at least 20 elements depthwise in 
the box {y direction ), and at least 40 in the x direction. All elements have a max- 
imum dimension of 0.06 dimensionless units. We have slabs which are at least four 
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Figure 13. Bout 


/ conditions used in the models. 
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elements thick. In cases where wa ara modelling the presence of discontinuities of 
soma kind (Sactlon 3.8), wa doubla this rasolution. 

Tha surfaca of tha Earth marks tha top boundary, and wa hava tha raasonabla 
conditions (Slaap, 1075; Rlchtar, 1073) that tha fluid is fraa to move horizontally, 
but has zaro vartlcal valoclty. Tha dapth 670 km. marks a bottom boundary which wa 
bagln by assuming marks an obstruction of soma sort to the motion of tha slab (a.g. 
Rlchtar, 1070). Wa thus set vartlcal valoclty to zero at this boundary. Horizontal 
velocity can ba sat to zaro or loft free. If horizontal valoclty is sat to zaro, tha boun- 
dary assumes tha character of an extrema viscosity jump. A fraa horizontal valoclty 
might ba batter for simulating a density discontinuity arising from chemical layering, 
whara lateral motions might ba occurring along tha deforming chemical boundary. Of 
course, vartlcal motions would occur also, In the dynamic situation. Wa ara obtaining 
Instantaneous solutions; wa cannot track tha dynamic deformation of a chemical 
boundary. Wa accept v v • 0 and v, fraa as an approximation to tha steady state at 
tha boundary In tha case of chemical layering. Leaving v M fraa Is a constant pres- 
sure boundary condition, with pressure equal to zero at the boundary. As we will soon 
see, the bottom boundary condition on v„ has very little Influence on our Important 
results— It Is tha condition on v y which Is important. 

The boundary conditions on the side of the box ara more arbitrary and more 
artificial. Since the Earth, or a part thereof, cannot be regarded as as Isolated box, 
It saems reasonable to allow fluid to enter and leavu our model box. Wa hava thus 
left tha horizontal valoclty fraa on the left hand side. In some models, wa will ba 
pushing the slab from the left. We have somewhat arbitrarily chosen to set v t equal 
to zero on the right hand side. It turns out that the boundary conditions on the right, 
while they may Influence some of tha details of the flow field, make essentially no 
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dlfference to the results for tho stresses In the slab, which arc our primary concsm. 
This can ba saan or. fig. 1 7, and figs. 01 7. 

Thara ara soma other assumptions wa are making. Wa are Ignoring Inertia and 
Coriolis forces, which ara much smaller in our problem than the viscous forces (Sieep, 
1076). Wa are also, obviously, ignoring the curvature of the Earth. As discussed by 
Richter (1073), this may cause some geometric distortion In the computed flow field 
compared to what might actually be happening in the Earth, but the basic dynamics 
will not be changed. 

There are many things we do not know about subducting slabs, and the subject 
can usually offer much room for debate. One thing everyone seems to agree on, how- 
ever, is that subducting slabs ara colder and denser than the surrounding mantle. 
Gravitational forces on the slab are likely to be an Important factor in the kinematics 
of the subduction process. How much colder one believes the slab to be than the 
mantle depends, of course, on one's thermal modal. McKenzie (1060) solved the heat 
conduction equation for a subducting slab In a mantle of constant temperature 7V 
His 50 km thick, 46° dipping slab moving at 10 cm/yr warms to roughly 0.8 7o at its 
coldest by the time It reaches a depth of 700 km. Hence, for T 9 ■ 1 300° , It Is over 
260° cooler than the surrounding mantle. Howard and Hager (1083) have performed 
a very Interesting calculation which refines McKenzie's model to Include the effect of 
the cooling of the mantle as well as the heating of the slab. The material parameters 
they use In their model are the same as McKenzie's. In their modal, the slab stays 
colder to a significantly greater depth than does McKenzie's slab, because of the 
buildup of a layer of cooler mantle on the sides of the slab. A 50 km slab in the 
Howard-Hager model warms up to 0.67 70 by 700 km; a 100 km slab warms up to 
0.46 7V Thus, for T 0 * 1300°, the slab may be about 700° cooler than the 
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surrounding mantle near the 670 km seismic discontinuity. 

Howard and Hager, like McKenzie, have not considered the adiabatic gradient in 
tha mantle, or the possible effects of shear heating and mantle phase transitions. 
Such effects have variously been considered in many studies: Mlnear and Toksoz 
(1970a,b);Toksoz et al (1971,1973); Hsui and Toksoz (1979); Turcotte and Schu- 
bert (1973); and Schubert et al (1976), among others. In several of these studies, 
the numerical calculations are so involved that it is sometimes difficult to isolate the 
effect of each Individual factor, and it is also difficult to compare the models to each 
other. The coolest slab appears to be that of Schubert et al (1 975), which at 660 
km or so is a maximum of 800° cooler than the surrounding mantle. A highly exoth- 
ermic olivine-spinel phase transition in the model may be partially responsible for this. 
We note that none of the models to which we have just referred treat the issue of 
the entrainment or the mantle, as Howard and Hager (1 983) do. 

Fig. 1 4(a) shows calculations of the density difference between slab and man- 
tle for the Howard-Hager Model with a 45° dipping slab moving at 6 cm/yr. The coef- 
ficient of thermal expansion is not well known, but is thought to be of order 10~ e 
(e.g. Sleep, 1976). We have somewhat arbitrarily adopted for a a value of 6 X 10~ s , 
and 3.6 g/cc for the density of the mantle at 0° C. As we can see from the figure, 
the density contrast varies across the slab and with depth. In most of the models, 
we have assumed the density difference between slab and mantle to be 0.07 g/cc 
throughout the slab. As we shall see in Section 3.7, this does not affect our results. 

Fig. 14(b) shows calculations of the viscosity ratio of the slab to the mantle. 

2£l ^ c R\ T 0 T l 
Vt 0 To* 


This ratio is of the form 
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Figure 1 4. Variation within the slab of (a) density contrast with respect to the man- 
tie and (b) logarithm of the ratio of slab viscosity to mantle viscosity. In 
both cases, the slab is dipping 46® and travelling 5 cm/yr. Calculations 
are from the model of Howard and Hager (1983). Vertical parallel lines 
Indicate boundaries of slab. Each curve is a perpendicular section through 
the slab, taken at the depth Indicated by the curve label. 
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whare T , la calculated temperature, T 0 la assumed Initial mantle temperature, E la 
the activation anargy for creep, and R la the universal gas constant. Wa have 
Ignored a term dependant on tha activation volume, since activation volume Is essen- 
tially unknown. Since tha precise composition of tha slab Is not known, wa can only 
guess at E. Wa use Ashby and Verrall's (1977) value for olivine, 6.2 X 10 4 J/mol. 
Tha figure shows that wa can expect viscosity In the slab to be everywhere essen- 
tially Infinite with respect to the mantle. 

Note, however, that the equation above calculates microscopic viscosity. 
Although we expect from the temperature that the slab Is absolutely rigid microscopi- 
cally, there Is evldance that Its macroscopic viscosity Is lower. Macroscopic and 
microscopic viscosities can differ If, for example, tha slab Is fractured. Melosh and 
Reef sky (1980) found that an effective viscosity of about 6 X 10** p Is required to 
explain the outer arc bulge and trench If these are formed by bending a viscous litho- 
sphere. The mantla Is believed to have a viscosity of roughly 10 M p from post- 
glacial rebound studies (Cathles, 1976; Peltier and Andrews, 1976). Thus our 
assumption of a slab viscosity ten times greater than the mantle viscosity Is probably 
reasonable. 

3.1.1 A Word About Units 

We solve a series of problems for a box which has depth h equal to 1 . The two 
fluids In the box simulating mantle and slab have viscosities of 1 and 10 respec- 
tively. We apply a downward body forca of 1 000 to the slab. The results we obtain 
ara not In any conventional unit. The reader will notice, however, that In all our fig- 
ures we report stress In bars and velocity In cm/yr. These results are correct when 
the problem is scaled such that h • 670 km, the viscosity of the mantle Is 10 8C p, 
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•nd a body fore* of 1 000 corresponds to a density contrast of 0.07 g/cc between 

slab and mantle. Stress scales as bpgh, and velocity scales as & ■ £*'' - , where g Is 

V 

acceleration due to gravity and 17 Is viscosity. It Is Important to note that stresses 
do not depend on the absolute value of the viscosity, whereas velocities do. To con- 
vert back to the "natural" dimensionless units of the problem, one needs only to 
divide stress In bars by 4.72, and no conversion Is required for the velocities. 

3.1.2 Presentation of Results 

We have calculated stressea and flow fields for a large number of models. For 
all these models, we present plots of stress magnitude versus depth In the slab (see 

3.2 below) In the main body of the text. In a number of cases, we also present spa- 
tial stress plots (showing stress orientation) and velocity fields; however, the large 
number of models would have made It Impractical to Include these In the text for 
every case. Since the spatial stress plots and the flow diagrams do contain Impor- 
tant information, we have placed the ones not Included in the text in Appendix D. 
The figures in Appendix D, as we explain there, are labelled so as to allow them to be 
easily identified with the corresponding figures In the text. 

3.2 A Sinking Vertical Slab 

We now examine the results of some simple calculations which yield much insight 
In explaining tha seismicity curves. The slab sinks under Its own weight, and Is not 
subjected to external push forces. Consider the model of a vertical slab, with 
parameters as discussed in the previous section. Fig. 1 6 shows the flow field. The 
slab Is moving at a velocity on the order of centimeters a year. Fig. 16 plots the 
axf»s of compression of the deviatoric stress at points throughout the model box (at 
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the center of aach flnita element). Tha ax as of tanslon era parpandicular to tha 
axas of compress'on, and of aqual magnltuda. Wa saa that tha slab Is In vertical 
(l.a., down-dip) tanslon on top, and down-dip compression on tha bottom. In fig. 1 7, 
wa piot tha average stress magnitude In tha slab versus depth. Wa define tha 
stress magnltuda as 


o = V 1/2 (t£+t*~+ Zrjy) 

where the T y are stress components. This Is actually the expression for tha second 
stress Invariant. At each depth, we calculate a for all points In the slab, take tha 
average, and plot on fig. 1 7. 

The calculated stress profile has tha following features: (1) a linear decrease 
down to w 300 km, (2) a minimum between 300 t.id 400 km, and (3) a resurgence at 
depth. In a rough sense, It seems to follow the same pattern as the curves of log# 
versus depth. The linear decrease In the upper part of the stress profile could 
explain the linear-log decrease of seismicity In the upper parts of subducting slabs 
quite nicely if the number of earthquakes depended exponentially on the stress, 
(e.g. Mogi, 1962a,b, and section 3.4). The position of tha stress minimum also sug- 
gests that the seismicity minima might be occurring at a depth dictated by the 670 
km length scale to a bottom barrier. We note that whether v t Is free or zero at the 
bottom boundary makes essentially no difference. The barrier can be a chemical one 


or a viscosity jump. 
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Flgure 1 6. Flow field for • vertical sinking slab extending to 070 km, subjected only 
to gravitational forces, as described in Section 3.2. 
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Figure 16. Deviatoric compression axes for the slab of figure 16 and section 3.2. 
Tension axes are of equal magnitude, and perpendicular in direction. 



Kbar 


116 - 


ORIGINAL PAGE IS 
OF POOR QUALITY 





Figure 17. Solid line shows average stress magnitude versus depth for the slab of 
section 3.2 and figure IS. Dashed line shows result for the same slab 
except that v x is left free on the bottom boundary instead of being set to 
zero. Dotted line shows result for the same slab when v m is left free on 
the right hand boundary Instead of being set to zero. 
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3.8 Viscosity Contrast at 670 km 

Suppose now that there is no barrier at 670 km. Fig 1 8 shows the stress profile 
for e slab extending to 670 km, but underlain by mantle fluid of viscosity 10 M p. (Our 
model box Is now twice as big in both dimensions). We do not see the minimum or the 
resurgence, but rather a smooth decay. The slab Is In tension, except at Its very tip. 

Fig. 1 0 shows what happens If we increase the viscosity of the lower mantle, 
creating a contrast at 670 km. As the* viscosity contrast Increases, the stress 
minimum, below which compression prevails, moves up in depth. A viscosity contrast 
of a factor of 6 produces a significant peak. A contrast factor of 26-60 produces a 
large peak, which is not much different In position or intensity from the case when 
the contrast is 1000— or essentially infinite, as in fig. 1 7. 

How viscous is the lower mantle in fact? At present, this question Is not com- 
pletely resolved, but there is some evidence that it could be more viscous than the 
upper mantle. A good discussion of the literature Is given in O'Connell (1877). As we 
have noted, the upper mantle appears, from analysis of post-glacial rebound data to 
have a viscosity of 1 0 22 p (There Is a possible low viscosity channel beir>w the litho- 
sphere; how low Its viscosity can be depends on Its thickness). Cathles (1876) has 
fit the data with a 76 km thick channel with viscosity o f order 10 zo p). The earlier 
literature (MacDonald, 1 863; McKenzie, 1 866) favored a large Increase of viscosity 
with depth in the mantle. McKenzie (1866) concluded that the lower mantle had a 
viscosity four orders of magnitude greater than the upper mantle. His conclusion was 
based on Interpreting the Earth's nonhydrostatic "fossil" bulge. However, It appears 
that his results were an artifact of his use of spherical harmonic coefficients. They 
were questioned by Goldreich and Toomre (1868), who placed an upper bound of 
10 s4 p on the viscosity of the lower mantle from the rate of polar wander obtained 
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Figure 1 8. Average stress magnitude versus depth for a vertical slab subjected only 
to gravitational forces and extending to 870 km, when there Is no barrier 
at 670 km. depth. 
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Figure 10. Average stress magnitude versus depth for a vertical sinking slab 
extending to 670 km when there is a viscosity contrast at 670 km depth. 
The number labelling each curve denotes the ratio of viscosity below to 
viscosity above the discontinuity. Viscosity above the discontinuity is 

1 0 88 p. 
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from paleomagnetlc studies. Cathles (1076) Interpreted tha Canadian post-glacial 
rebound data to Indicate a lower mantle viscosity of 10 M p the same as the upper 
mantle. Walcott (1073) believed the data to indicate a lower mantle viscosity of at 
least 10* 9 p. As O'Connell (1077) points out, Walcott's analysis considered the resi- 
dual gravity anomalies In Canada associated with unrecovared rebound. Cathles 
(1076) believed the data to be consistent with either his or Walcott's model, the 
choice depending In large measure on the significance attached to the residual grav- 
ity anomalies. More recently, Yuen et al. (1082) estimated the viscosity of the lower 
mantle to be larger than that of the upper mantle by analysing observed secular 
motions of the Earth's rotation axis. They found it to be larger than the viscosity of 
the upper mantle, at most by a facta* of 4. Hager (1083), by considering geoid 
anomalies, has found that the contrast factor must be as high as 30 In subductlon 
tones. Thus we see that, overall, tha available data are consistent with enough of a 
viscosity barrier at 870 km to produce stress patterns in the slab matching observed 
seismicity profiles. 

8.4 Relationship of Seismicity to Stress Levels 

Before exploring more models, we pause to consider the important question of 
the relationship of seismicity to stress levels. Qualitatively, our results make sense If 
we assume only that a higher level of stress leads to a larger number of events. 
They make even more sense if the number of events depends exponentially on the 
stress. The linear-log decease of seismicity In the shallow portion of slabs, as well 
as the relative numerical levels of deep and shallow seismicity, are nicely explained. 
All this assumes that the dependence of seismicity on stress does not change drasti- 
cally with depth. 
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Let us furthsr examine the Mm that N a e* w , whare k Is soma constant, ano cr 
danotas the stress magnitude. It has long been known (Ishlmoto and llda, 1 939) that 
earthquakes follow a frequency-magnitude relation of the form relation of the form 

k>oV= a-bli 

where a and 6 are positive constants known universally In the seismological litera- 
ture os the "a value" and the ”b value” (Richter, 1968; Bath, 1981). We may 
rewrite this as 


togN=b(M nMX -M) 

where the maximum magnitude of earthquake Is equal to £ — provided, of 

course, that the linear relationship holds throughout the magnitude range. If b is not 
a function of stress, we have that N a e*“*. If M m a c, then M a **'. There Is no 
compelling reason to assume that Mm a a. We are merely trying to see the condi- 
tions under which the empirical distribution of earthquake sizes might lead to the type 
of exponential dependence on stress which has been observed experimentaly, and 
which would provide a link between our slab stress calculations and observed seismi- 
city curves. 

The experiments where N a a kc has been observed have been acoustic emis- 
sion studies on the brittle failure of rock samples. In evaluating such experiments in 
connection with our problem, we are faced with the omnipresent difficulty of scaling. 
How relevant rock mechanics experiments on small samples are to the real Earth is a 
long standing unresolved question (e.g. Ilo, 1 982). We are also faced with the diffi- 
culty of applying results from brittle fracture experiments to deep focus earthquakes 
whose mechanisms, albeit consistent with shear dislocations, are not known to result 
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from brittle falure. However, the reeulte of the experiments are worth reviewing. 
This b because, first of ell, these results ere al we now have at our disposal. 
Secondly, we must bear In mind that brittle failure, although not a certain mechanism 
for earthquakes at depth, Is still a possible one, particularly If the stab retains signifi- 
cant amounts of pore fluid. This is partlculariy true for the relatively shallow depths 
(above 260 km) where we are interested In finding an explanation for the linear-log 
decrease in number of events with depth. 

Mogi (1862a,b) performed experiments establishing that microfractures In rock 
samples obey a frequency-magnitude relation similar to that of earthquakes. He 
found in his experiments that the number of events was proportional to the exponen- 
tial of the applied stress. This dependence was subsequently used in certain seismic 
hazard studies (Hagiwara, 1974). Scholz (1968) conducted more rock mechanics 
experiments, and found the b value to decrease with stress. This complicates the 
exponential dependence of number of events on stress. However, this stress depen- 
dence of the b vaiuo has subsequently been disputed. Mogl (1981) presents data 
showing that under a constant load, the 6 value of the microfiactures decreases 
before the main rupture. Ohnaka and Mogi (1981) argue that "the relative increase 
In the number of emission events with larger amplitude [ l.e., the decrease in the b 
value] is not due directly to the increase In the stress level Itself. The effect can be 
explained reasonably If the stress drop and/or the source dimension (crack size) of 
emission events become larger in statistical terms as the rock approaches failure". 
Other data (Kusunose at al., 1980; Sano at al., 1982) have also displayed the effect 
of a relatively constant b value until the sample comes dose to failure. The b values 
of Kusunosa et al. (1980) are constant until the stress exceeds 85% of the fracture 
strength. Sano et al. (1982) believe that their b values show evidence of 
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decreasing •trass, but this Is, again, for a deformation regime closa to tha macros* 
copic failura of tha rock. Thara Is thus no convincing experimental evidence that tha 
b value changes significantly with stress. 

If the 6 value did change with stress, then we might expect from our calcula- 
tions that it would change with depth In the Earth. Gutenberg and Richter (1964) 
reported b values for shallow, Intermediate, and deep earthquakes respectively as 
0.9 ± 0.02, 1.2 i 0.2, and 1 .2 ± 0.2. Fig. 20 shows frequency-magnitude plots for 
the catalogs we have used in this study for shallow, Intermediate, and deep earth- 
quakes. We focus our attention on the magnitude range between 6 and 6, where the 
linearity of the frequency-magnitude relation is most apparent For the smaller 
events the relationship breaks down, because of inadequate detection capability. 
For larger events, It breaks down again, both because of the saturation of the m* 
scale (Kanamorl, 1978) and possibly because of Inadequate sampling of larger 
events over the sixteen year time period. There Is no evidence In fig. 20 that 6 
value changes with depth. Plotting a figure like fig. 20 for 1 00 km depth Intervals 
also falls to provide such evidence. Chouhan and Srlvastava (1970) report a depth 
dependence of b value based on an analysis of Gutenberg and Richter's (1964) data 
In 100 km depth Intervals. Their results indicate a roughly constant b value of 0.66 
for events above 400 km in depth, a value of nearly 0.8 below 600 km, and a value 
of almost 1.1 from 400-600 km. We do not believe these results are significant 
given the small number of events In each sample. The observation by Abe (1981) of 
a 6 value Increasing with depth for the large events In the data set of Abe ar.d 
Kanamorl (1979) suffers from the same limitation, as we can see by examining fig. 
21. The b value for the deep events Is not well constrained. We note, finally, that 
Kagan and Knopoff (1980), In a rigorous statistical study of the NOAA catalog, find no 



Figure 20. Logarithm of the number of earthquakes versus one-second body wave 
magnitude m*, In Intervals of 0.1 magnitude unit A is for shallow (less 
than 70 km) events; B Is for Intermediate (betwoen 70 and 300 ;m) 
events; and C Is for deep (greater than 300 km) events. Data sources 
are the NOAA and PDE catalogs, 1064-1080. 
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Figure 21 . Logarithm of numbar of earthquakes versus long partod body wava magni- 
tude mg for large events of the twentieth century compiled by Abe and 
Kanamorl (1070). Figure 21(a) is for deap events; figure 21(b) is for 
intermediate events. 
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signlflcant differences in the 6 values of shallow, intermediate, and deep events. 

3.6 Finite, Dipping Slabs 

All slabs do not extend to 670 km, and we might raise the question of whether 
finite slabs in our model can give us acceptable stress profiles. Fig. 22 shows that 
this is in fact the case. A vertical slab extending to less than 300 km, sinking under 
Its own weight, has a smooth decay of o with depth. One extending deeper develops 
a peak, as the bottom boundary begins to be felt. We note that the shallow slab of 
fig. 22 is in tension, while deeper ones are In tension on top and compression on the 
bottom below the stress minimum. 

What happens when the slab is dipping? Fig. 23 shows the results for a 60° 
slab, and fig. 24 for a 45° slab. AH parameters in the models are the same as for the 
vertical slab, except that the slab is now dipping. The overall stress levels are lower 
than in the vertical slab because the body force is the same, but is not directed 
down-dip. Overall, the dipping slabs show the same qualitative features as the verti- 
cal one. There is more difference between the v M free and the v x =0 case when the 
slab extends to 670 km, but this is still not a difference resolvable in the seismicity 
data. 

In general, the shallow decay, the minimum, and the resurgence are there. There 
le some complexity caused by the dip, because of sagging effects. The shallow 
decay is less smooth, the deep peak is not as pronounced, and the minimum is some- 
what broader In the 46° case than in the 60° case, and this is also true of the 60° 
slab with respect to the vertical one. Of course, togAf versus depth for the various 
regions is not always smooth in the shallow portion, and does not always show a pro- 
nounced deep peak. We cannot, however, see any pattern of the complexity of the 



Figure 22. Stress profile for s vertical slab subjected only to gravitational forces, 
extending to (a) 270 km, (b) 400 km, and (c) 640 km. 
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Figure 23. Stress profile for a 00* dipping slab extending to (a) 270, (b) 402, (c) 
640, and (d) 670 km. Dashed curve in (d) is for v, free on the bottom 
boundary; solid curve is for v, =0 on the bottom boundary. Figure 23(d) is 
the only case where the two curves do not effectively coincide. The slab 
Is subjected only to gravitational forces. 
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Figure 24. Precisely as figure 23, but for a 46 s dipping slab. 
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salsmic profiles varying with the dip. It is. In fact, almost Impossible to Identify any 
trends In the seismicity profile which we might ascribe solely to the dip. 

The stress orientation for a 45° slab Is shown In fig 25. Unlike when the slab Is 
vertical, the stresses In the depths above 300 km are not precisely down-dip tensile. 
There are both tensile and compressive components down the dip. The state of 
stress In the lower portions of a deeply extending slab, however, Is clearly down-dip 
compressive. This is consistent with the famous conclusions of Isacks and Molnar 
(1071), and also with our own analyses of earthquake stress axes obtained by 
moment tensor inversion (Section 2.2.2, Appendix A). 

3.0 Slabs Pushed from the Side 

The models we have been considering so far have ail been subjected to a body 
force only. While gravitational forces on the slab are likely to constitute an Important 
driving mechanism for plate tectonics (e.g., McKenzie, I960), ridge-push is also a 
possible force. Fig 26 shows a series of stress profiles for a slab pushed from the 
side so as to produce stress levels in the slab comparable to those produced by body 
forces in our previous models. We see that we can match the smooth decay in the 
shallow portion, but we c**i not produce a deep peak. As shown in fig 27, stress is 
down-dip compressive throughout the slab, in general, a model with push force only 
Is not adequate to satisfy the observations. We have no trouble, however. If we add 
a push force to a body force, as shown in fig 28. As we can see In fig. 29, the state 
of stress in the shallow portion of the slab Is less down-dip tensile when we add a 


push force. 
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Figure 26. Compression axes for the 46* slab of figure 24, extending to (a) 270 km, 
(b) 400 km, (c) 640 km and (d) 670 km. 
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Figure 26. Stress profiles for a 46* slab subjected only to a push force from the left 
side. Different curves are for slab extending to different depths: 270, 
400, and 670 km. 
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Figure 27. Compression axes for the slab of figure 26 extending to 670 km. We 
recall that the siib is dipping 46 * and is not subjected : r Qvitationai 
forces, but only to a push from the left side. 
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Figure 28. Stress profiles for a 45* slab axtendir .o various depths (270, 400, 
540, and 870 km). The slab is subjects ,.oth to the push force of figure 
26 and to a gravitational body force as i, figure 24. 
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29. Cumvress^n a xa- for the slao of figure 28 extending iO 670 km. We 
recall that this slab is dipping 45° and is subjected both to a push from 
the side and to gravitational body forces. 
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3.7 A Slab Mora Consistent with Thermal Models 

As we have noted, the assumption of a uniform density throughout the slab is 
somewhat oversimplfied. If we look at fig 14, the calculation of density differences 
between the slab and mantle for Howard and Hager's (1983) model, we see that bp 
varies both across the slab and with depth. Fig 30 shows a stress profile for a slab 
with body force decreasing smoothly with depth, and one for a slab with a laminated 
body force. Clearly, we may regard constant density throughout the slab as a rea- 
sonable approximation for our purposes. 

3.8 Discontinuities and Phase Transitions Above 070 km 

All our slabs thus far have been descending in a uniform mantle with the barrier 
at 670 km being the only discontinuity. There are known seismic discontinuities 
between the Moho and 670 km, however. One is the Lehmann discontinuity (Lehmann, 
1 961 ; Hales et al., 1 976) In Anderson's (1979a,b) compositional model of the mantle, 
this discontinuity represents a chemical boundary between peridotite and eclogite. 

Fig 31 shows a stress profile for a 45° slab when there is a barrier to vertical 
flow at 200 km. Horizontal flow Is allowed, as discussed in Section 3.1 for a chemical 
boundary. The barrier does not extend into the slab. Hence, we leave v a free and 
set v y equal to zero at the internal surface representing the barrier, and that surface 
does not include the 1 1 element nodes representing the slab. We see that a chemi- 
cal discontinuity in the mantle would be expected to produce a peak in stress. Thus, 
within the framework of our assumption that the seismicity profile follows the stress 
profile, there is no evidence for a barrier of this sort above 670 km. 

Another discontinuity is the bitter known one at 400 km (e.g., Anderson, 1967). 
This second discontinuity has been widely thought to represent a phase transition of 
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Figure 30. Stress profiles for slabs more consistent with the results of the thermal 
model shown in figure 14. The solid line is for a slab dipping 45° and 
extending to 670 km, with a dens it) structure such that the upper surface 
of the slab is densest, with a gradient down to the lower surface which is 
the least dense. The dashed line is for a similar slab, but with density 
decreasing with depth to half its shallow value. 
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upper mantle divine Into the spinel structure (Meijerlng and Rooymans, 1968; Ring- 
wood and Major, 1970). The proposed olh/lne-spinel change has been thought to 
have significance In plate tectonic processes, In that Its positive Clapeyron slope 
Implies an elevation of the phase boundary in the cold slab, this possibly providing a 
significant additional downward body force (Turcotte and Schubert, 1971). Schubert 
et al (1975), have calculated this body force in their model descending slab. They 
estimate a quadrupling of the body force above the 400 km discontinuity. Fig 32 
shows the stress profile for a 46° slab where this effect is included. The influence 
of this increase in body force on the stress profile is very significant, and no evi- 
dence for such behavior Is seen In the seismicity profiles. 

Thus we affirm once more that the distribution of seismicity with depth in the 
world appears, to first order, to be dictated by the presence of a barrier at 670 km. 
Other proposed chemical discontinuities and phase transitions do not appear to be 
necessary to account for the observations, and in fact show evidence of being 
Inconsistent. Phase transitions where the phase boundary is not significantly 
elevated in the slab may be consistent. It Is interesting to note that the phase tran- 
sition of clinopyroxene to the garnet structure, favored by Anderson (1979a) for 
explaining the 400 km discontinuity, may have a Clapeyron slope lower than that for 
olivine -* spinel, although the value is not well constrained at present (Akaogl and 
Akimoto, 1977). It should be borne in mind, however, that for any phase transition, 
kinetic factors might prevent a significant change in elevation of the phase boundary 
In the slab (Hager, personal communication, 1 983). 
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Figure 31. Stress profile for a 46° slab extending to 670 km, but with a barrier to 
vertical flow at 200 km arising from a hypothetical chemical discontinuity. 
The barrier does not extend into the slab, which in this scenario is 
assumed to have pierced the discontinuity. 
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Figure 32. Stress profile for a 46° slab with an elevated olivine-spinel phase boun- 
dary as calculated by Schubert et al. (1075). A is for a phase boundary 
elevated by a 100 km, as In that paper, and B is for a phase boundary 
elevated by 50 km. 
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8.9 Viscocity Increasing with Depth in the Mantle 

If we do not Impose a step function Increase in viscosity at 670 km, but rather 
allow viscosity to Increase smoothly with depth in the mantle, do we get the same 
results? Fig. 33 shows the results of a model where viscosity in the upper mantle 
Increases linearly with depth such that the viscosity is 26 times greater at 670 km 
than It Is immediately below the lithosphere 'where It is 1 , that Is, 1 0* z p in our scal- 
ing). The lower mantle, below 670 km, has a viscosity of 25. The viscosity cf the 
slab Is kept at a value of 10 as in previous moaels. The stress profile H smoothed out 
somewhat; the deep peak is not as intense. However, the profile still shows some of 
the first order features observed in the seismicity. 

This Is less true of the profile In fig 34. We arrive at this profile by using the 
model of fig 33, except that we now allow viscosity to Increase with depth in the 
slab as well as in the mantle. The slab is now ten tlmos more viscous than the mantle 
at every depth. There is no reason to presume. If viscosity Increases with depth in 
the mantle, that it should not also do so in the slab. Thus we might say that, based 
on the poorer match of fig. 34, an increasing-viscosity model is not as successful in 
explaining the variation of seismicity with depth as Is a sharp increase at 670 km. 

3.10 The Asthenosphere, Mesosphere, and Double Benloff Zones 

We explore here the effect of a low viscosity zone in the mantle just below the 
lithosphere (Cathles, 1 876). We refer to this 2 one, to which we assign a viscosity of 
0.01 .slatlve to the rest of the mantle, as the "asthenosphere", following established 
convention, and the more viscous mantle below it as the "mesosphere ", following 
Sleep (1879). The sagging stresses In the slab as it meets the asthenosphere- 
mesosphere contact have been proposed (Sleep. 1978) as a possible explanation of 
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Figure 33. Stress Profile for a 45° dipping slab subjected to gravitational forces 
only, when viscosity in the mantle increases linearly with depth from a 
value of 1 (10 82 p) immediately below the lithosphere to one of 25 at 670 
km. The viscosity of the slab is 1 0 at all depths. 
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Figure 34. The same model as that in fig. 33, except that viscosity increases with 
depth in the slab as well as in the mantle, so that the viscosity of the slab 
is ten times that of the mantle at every depth. 
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double Benioff zones. We will see below that models of the kind considered by Sleep 
also yield stress profiles consistent with observed seismicity. 

The literature on double Benioff zones has been reviewed by Fujita and Kanamori 
(1081). A double Benioff zone is a structure of seismicity involving two parallel 
planar earthquake belts, separated by about 30-40 km, and existing between 
approx'mately 65 and 1 85 km depth. The upper zone is characterised by down-dip 
(or, as Fujita and Kanamori call it, "in-plate") compression, while the lower zone is 
tensile. This stress orientation is opposite to what one would expect for a bending 
plate, and in fact "unbending" of the .4,'ate was an early proposed explanation for the 
observations (Engdahl and Scholl, *,977). Fujita and Kanamori cast some doubt on 
the unbending hypothesis, because they argue that if this were the explanation, one 
might expect double Benioff zones to be a global feature. In fact the only true dou- 
ble Benioff zones known are in Japan (Tsumura, 1973; Umino and Hasegawa, 1975, 
1982; Hasegawa et al., 1978) and the Kuriles (Veith, 1974, 1977). Engdahl and 
Scholz (1977) reported a double Benioff zone for the Aleutians. However, the 
existence of this zone has been called Into question by Topper (1978, cited by 
Fujita and Kanamori), who believes that what one is actually seeing in the Aleutians is 
a tear in the Benioff zone, whioh sh ms up as an extra plane when all everts are pro- 
jected onto a vertical cross section. Reyners end Coles (1982), who continue to 
treat the Aleutian data as a double Benioff zone, report that the stress orientation in 
the upper and lower planes is opposite to what one might expect from unbending, 
that is, opposite to what is seen in Japan. 

Sleep's (1979) models have a lithosphere of 2X10 zs p viscosity, an astheno- 
sphere of 2X1 O zo p, an accretionary wedge of '.cosity varying between 2X1 0 a: and 
2X10 as p. and a mesosphere of viscosity varying between 1 0 2 1 and 2X10 !2 p. The 
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accretionary wedge does not have much effect on Sleep's results, but the meso- 
sphere does. If mesospheric viscosity Is less than about 2X1 0 zz p, the stress in the 
slab is down-dip tensile. If it is higher than 5X10 ez p, the stress is down-dip 
compressive. Between these two viscosities, a double Benioff zone develops. 

Fig. 36 shows a calculation very similar to Sleep's for a mesospheric viscosity of 
2X10 zz p. (Note that Sleep's boundary conditions are somewhat different from those 
we have used in our other models. He imposes a velocity of 5.7 cm/yr on the litho- 
sphere on the left side. We have applied similar boundary conditions). The slab 
extends to 230 km depth and dips 00® — Sleep was trying specifically to explain the 
situation then thought to prevail in the Aleutians. The top of the mesosphere is at 
160 km. Examining fig. 36(b), we see that below about 100 km, the upper plane of 
the slab is in compression, and the lower plane Is (somewhat weakly) in tension. Fig 
35(a) plots the stress magnitude as in our other calculations. It behaves like our 
other shallow profiles, showing a fairly smooth decay. If we let the slab extend to 
670 km, '••ith our customary barrier at that depth (fig. 36), the stress profile is not 
qualitatively different from that in a uniform mantle. The minimum is just somewhat 
broader. At depths below about 400 km, the stress is down-dip compressive 
throughout the slab, while at depths from about 120-400 km, the stress is as in a 
double Benioff zone. This is deeper than the range associated with double zones. 

3,11 Non-Unear Rheology 

All calculations so far in this study have assumed tliat the mantle and slab 
behave as Newtonian fluids. Other than computational convenience, there is no 
a priori reason to assume this. While some proposed deformation mechanisms for 
the mantle, such as point defect dlf'usion (Nebarro, 1948; Herring, 1950; Raj and 



Figure 35. (a) Stress profile for a 60° slab extending to 230 km, with an 
asthenosphere-mesosphere contact at 160 km depth. The viscosity of 
the asthenosphere is 0.02 X 10 2t p, and of the mesosphere 2 X 10 Z2 p. 
Figure 35(b) shows compression axes for ’'his slab. 
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Figure 36. As Figure 35, but for a slab extending to 670 km. (a) plots the stress 
profile, and lb) plots the compression axes. 
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Ashby, 1071) or superplastic creep (Twlss, 1076) have linear constitutive relations, 
others do not. Diffusion controlled dislocation climb, known as "dislocation creep", 
has a strain rate proportional to some power n of the stress, n can be between 3 
and 10, tending to the higher end of the range with increasing temperature and 
decreasing stress (Stocker and Ashby, 1073). 

We have calculated stress profiles for models with a nonlinear rheology. As n is 
not known, we have chosen a value of 3 for computational reasons. The results for a 
46° dab extending to 670 km and subjected only to gravitational forces are shown 
in fig 37. They are very similar qualitatively to the linear results, and wo find no 
cause to suspect that nonlinear mantle rheology will alter or Invalidate our conclu- 
sions. We note that the calculation assumes a slab viscosity 100,000 times greater 
than the mantle viscosity of 1 0 88 p at unit stress (4.72 bars, as discussed in Section 
3.1.1). The contrast In effective viscosity, however. Is far lower— viscosity is 
stress dependent. We can see this In fig. 38, which snows contours of effective 
viscosity, where a viscosity of 1 represents the starting mantle viscosity. The siab 
has an effective viscosity of order 10, while the mantle varies from order 0.01 to 1. 
The nonlinear rheology tends to create a low viscosity zone where the slab Is bend- 
ing, and on the underside near the bottom where the slab approaches the barrier. We 
can see hints of a tensional lower surface and a compressive upper surface in the 
dab (fig. 037S), as we observed in our previous linear models containing an astheno- 
sphore (Section 3.10). 

We note that the choice of the factor 10 9 was made in order to obtain reason- 
able flow velocities cf order 1-10 cm/yr. In the linear case, this is achieved by 
choosing t) SJ jW = 1 0rjyjwriE- * nd contrasts of 1C 9 would lead to very small veloci- 
ties. in the nonlinear case (n = 3), choosing rfsUB - 1 0 V° mantle Itads to 
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Figure 37. Results for a nonlinear rheology (n ■ 3). Slab dips 45° and Is subjected 
only to gravitational forces. Viscosity of slab at unit stress is 1 0 s times 
that of mantle. 
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Flgura 38 Contours of cffactlvs viscosity for the slab of fig. 37. A viscosity of 1 
ropressnts ths starting mantis viscosity. Contours ara drawn at affsctlve 
vlscosltlss of 0.06 ,0.1 ,0.3, 1 .0, 3.0, 10.0, and 30.0. 
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velocities on the order of tens or hundreds of me tars e year. 

3.1 2 The Nature of the 670 km Barrier 

Whether the 670 km harrier is a chemical discontinuity or a viscosity jump is not 
resolvable by our data. The stress patterns— both orientation and magnitude— are 
explainable , to first order, either way. The flow patterns, however, may provide 
some constraint. If we look at fig. 38, the flow field for a slab running into an 
impenetrable barrier, we see that the flow lines near the bottom of the slab do not 
parallel the dip of the slab, but rather flatten out. The slab looks as If it is about to 
distort into a sigmoidal shape. This kind of shape is not observed. Looking at fig. 40, 
the flow lines for the slab when there is a viscosity jump of an order of magnitude at 
870 km, we see that the flow roughly parallels the slab, and no distortion is implied, 
particularly when the slab Is allowed to extend into the lower mantle. This result is 
related to that obtained by Hager and O'Connell (1 978). They calculated global man- 
tle flow fields driven by observed surface plate motions. They were able to predict 
subduction zone dip angles successfully when they allowed flow to extend into the 
lower mantle, but not when they confined flow to the upper mantle. 

There are other arguments favoring the viscosity contrast as well. Some of the 
most powerful are the ones given by Hager (1983). He finds that geoid anomalies 
over subduction zones are positive, contrary to what one might expect in a model 
Earth of uniform viscosity; a factor of 30 or more Increase In viscosity with depth is 
required to account for this observation. Moreover, Hager finds that the density con- 
trasts in the seismically active parts of Benioff zones are insufficient to account for 
the magnitude of the observed geoid anomalies. Allowing slabs to penetrate aseismi- 
cally into the ower mantle can account for this discrepancy. Another way to account 



188 - 


for It, under the hypothesis of a chemically layered mantle. Is to have 360 km of dead 
slab at the base of the upper mantle. Hager argues that this would require a sub- 
stantial deflection of the 670 km discontinuity, which has not been observed. 

The finer 't atures of the seismicity patterns In some of our regions may be 
better explained by a viscosity contrast of a factor of 1 0 to 30 than by a very large 
viscosity contrast or a chemical discontinuity, if we refer back to fig. 19, and then 
to the plots in Appendix C, we see that many of the deep regions have deep peaks 
whose size relative to shallow seismicity levels might imply a "soft" boundary at 870 
km. The Marianas, for Instance, have a deeply positioned seismicity minimum, and a 
sr.mll deep peek. This pattern is better matched by a vertical slab sinking into a 
viscosity contest of an order of magnitude , or even less, than by one sinking into a 
hard boundary. The seismicity pattern ooserved In Mindanao is also suggestive of a 
relatively soft boundary, while that seen in Tonga Is more suggestive of a hard one. 
Perhaps we are seeing the effects of lateral variations in mantle viscosity. 

So can we conclude that the boundary at 670 km is a viscosity contrast rather 
than a chemical discontinuity? Not necessarily; consider fig. 41, which plots the 
stress profile for the slab of fig. 40(b). The upper mantle stress pattern Is adequate 
to explain observed seismicity, but the stress In the slab Is still high at depths below 
670 km. Why, then, do earthquakes not occur there? One can adopt the somewhat 
ad hoc explanation (e.g., Wortel, 1982) that w 700 km simply happens to be the 
depth at which a critical temperature Is reached and the slab loses Its mechanical 
integrity, becoming aseismlc. Such a:i explanation is not satisfying, however; it 
seems like too much of a coincidence that earthquakes should stop, after a very high 
level of activity, at a depth where there Is also a possibly very sharp seismic discon- 
tinuity (Richter, 1979). Thus, the question is not resolved. 



-188 


The sharpness of the discontinuity (Whitcomb and Anderson, 1968) has baen 
used by Anderson (1978) as an argument against Its being caused by a phase tran- 
sition. Bell (Geological Sciences Seminar, Caltech, 1 983), however, has disputed this 
argument, stating that the transition of upper mantle phases to the perovsklte struc- 
ture may indeed be sharp enough to explain the discontinuity. Thus, a phase transi- 
tion is stm a possibility, so we examine the possible effects of one kind of transition 
which has the potential to help us explain the seismicity patterns. If a phase transi- 
tion at 670 km has a negative Clepeyron slope, the phase boundary may be 
depressed In the slab. This may then cause an upward body force on the slab, pro- 
ducing compressive stress at depth, and possibly a stress profile with a minimum and 
a deep peak. It was once believed that upper mantle olivine In the spinel structure 
might disproportionate into mixed oxides (Birch, 1962), a transition which might 
Involve a negative Clapeyron slope (Ahrens and Syono, 1967). Since the discovery 
that olivine and pyroxene transform to the perovsklte structure at high pressure (Liu, 
1976,1976), there is less reason to suspsct a negative Clapeyron slope. However, 
It is still possible, as the slope for the perovsklte transition Is not known. 

Fig. 42 shows a curve for the spinel to mixed-oxides transition as calculated by 
Schubert et al. (1976), with a density contrast of 0.4 g/cc, and a « 30 km depres- 
sion of the boundary In the slab. We see that a model of this sort is capable of ade- 
quately matching the seismicity profiles, without an Increase In mantle viscosity or a 
chemical barrier at 670 km. However, as we have said, we have no reason to sup- 
pose a negative Clapeyron slope. A viscosity change or a chemical discontinuity still 
constitute more straightforward explanations of the phenomenr at hand. 



180 - 


Flguta 38. Flow field for a 45* slab sinking under Its own weight whan there Is a bar- 
rier at 670 km. The stress field for this slab Is shown in figs. 24 and 25. 
The flow lines appear to flatten out at the boundary, suggesting that the 
slab may be about to distort Into a sigmoidal shape. 
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Figure 40. Flow fields for a 45° slab sinking under its own weight when there is a 
viscosity contrast of an order of magnitude at 670 km. In (a), the slab 
extends to 670 km depth. In (b), it penetrates into the lower mantle, to a 
depth of 1 000 km. The flow lines are more aligned with the dip of the slab 
than they are in fig. 39. 
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Flgure 41. Stress profile for the slab of fig. 40(b), a 46° slab sinking under its own 
weight Into the lower mantle to a depth of 1000 km, when there Is a 
viscosity contrast of an order of magnitude at 670 km. The stress in the 
upper mantle portion of the slab Is consistent with observed seismicity, 
but we note that the stress In the lower mantle portion is still high. 
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Figure 42. Stress profiles for a 45° slab sinking under its own weight when there is 
a phase transition at 670 km with a negative Clapeyron slope, as calcu- 
lated by Schubert et al. (1975). 






4. Summary and Conclusions 


(1) . The distribution of earthquakes with depth in the world has the following 
features: (a) an exponential decrease from shallow depths down to a 250 to 300 
km, (b) a minimum leve! near 260 to 300 km, and (c) a deep peak below 300 km. 
Many shallow subducting slabs show only feature (a). Oeeper extending regions 
tend to show (a), (b), and (c), with the deep peak varying in position and intensity. 

(2) . A survey of events analysed by moment tensor inversion has confirmed 
some earlier Ideas concerning the state of stress in the slab. Deep earthquakes 
(below 300 to 400 km) tend to have compression axes aligned with the dip of the 
slab. This appears to be a global tendency. The behavior of intermediate earth- 
quakes is less clear, and more region-dependent. Both shallow-extending slabs and 
deep-extending slabs other than Tonga have intermediate earthquakes which are, in 
orientation, closer to down-dip tension than they are to anything else, but whose 
tension axes are not as well aligned with slab geometry as are the compression axes 
of deep events. The Tonga region shows some tendency toward down-dip compres- 
sion at intermediate depths. In general, however, we do not agree with earlier con- 
clusions ((sacks and Molnar, 1971) that deep-extending slabs in general are in 
down-dip compression at all depths. 

(3) . Simple viscous fluid models of subduction can explain observations (1) and 
(2) very well if there is a barrier of some sort at 670 km depth. A wide variety of 
models have stress magnitudes in the slab which display the following features as a 
function of depth: (a) a roughly linear decrease from shallow depths to about 250- 
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300 km, (b) • minimum near 260*300 km, and (c) a daap paak whose position and 
Intanalty depend on the dapth of penetration if tha atab. If the 670 km dapth marks 
a viscosity contrast, tha position and Intanalty of tha daap paak ara also affactad by 
tha magnltuda of tha viscosity Increase: "softer" boundaries produce deeper and 
smaller peaks. Curves of stress magnltuda versus dapth look vary much Ilka curves 
of log N versus dapth, where N Is tha number of earthquakes. Tha minimum at * 300 
km seams to be dictated by tha 670 km length scale. Tha linear decrease In log N 
with depth down to * 300 km Is understandable If the number of events depends 
exponentially on tha stress, for which there la some experimental evidence. Slab 
models with a barrier at 670 km yield down-dip compression below 300 to 400 km, In 
concordance with observation (2). At Intermediate depths, dipping slabs <n tha 
models have a state of stress which Is nelthar down-dip tensile nor down-dip 
eompreeslve, although It Is closer to the former. This, too, agrees with (2). 

(4). The observations are explainable If the slab sinks under Its own weight and 
Is not subjected to push forces. They are alw explainable If push forces exist In 
conjunction with gravitational forces. However, a slab subjected to push forces 
alone does not develop a deep peak In stress, and thus Is not adequate to explain 
the observations. 

(6). Chemical discontinuities above 670 km, or phasa transitions with a phase 
boundary elevated In the slab relative to the mantle, produce peaks In stress which 
do not appear to be mirrored In the seismicity. Phase transitions whose boundaries 
are not elevated In the slab may be consistent with the observations. 

(6). The date are consistent with s uniform viscosity mantle above the Impoi - 
tent barrier at 6 "0 km. Inclusion of a low vlacoalty asthenosphere below the litho- 
sphere does not destroy the match betwaen calculated stress profiles and observed 



seismicity, and oan also help explain the atraaa orientations aaaoclatad with double 
BanlofT zones (Sleep, 1070). 

(7) . The results for a nonlinear rheology (n ■ 3) are qualitatively similar to the 
linear results. 

(8) . Observations (1 ) and (2) ara equally well explained If the barrier at 870 km 
la a chemical discontinuity or a viscosity contrast where viscosity Increases by an 
order of magnitude or store. A viscosity contrast yields flow fields In the models 
which are more consistent with observed slab shapes. However, If we allow the slab 
to penetrate Into the lower mantle, we find that stresses In the slab below 670 km 
are as high as they are at upper mantle depths. We are thus faced with the problem 
of explaining why these stresses do not produce earthquakes. 
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Appendix A: Mechanisms of Larga Intarmadiats and 
Daap Earthquakes, 1978-1081 


Ou 

We present here the results of surface wave inversions for the focal mechanism 
of large intermediate and deep earthquakes from 1978-1981. Table A1 lists infor- 
mation from the NOAA catalog about these events, in deciding which events were 
'large", we were guided by the magnitudes. We studied earthquakes with m b ^6, and 
which offered some other evidence of being large, for example, a large listing under 
Af oMar . Magnitudes, however, are not necessarily reliable indicators of the size of an 
earthquake (Kanamori, 1978). Thus, some of the events we studied did not in fact 
turn out to have moments which would place them among the very largest events for 
the period. 

Our data consist of vertical component Rayleigh waves recorded on IDA (inter- 
national D enjoyment of Accelerographs) instruments. Generally, we use R2 in the 
inversions. Often, we also use R1 , and sometimes R3. The technique is described by 
Kanamori and Given (1981,1982). What follows in this paragraph is taken almost 
directly from Section 2 of the 1981 paper. The parameter vector is 
M-lMgy,M n -M„,M V y+M m ,M y9 ,M aa '} T . Once M is determined, the moment tensor 
matrix 



Mgy 



Myy 

My, 


My, 

M„ 


can be diagonaiised. The trace of the matrix is zero by assumption in the inversion. 
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The eigenvectors daflna tha principal stress axes. If tha Intermediate axis Is zero, 
the moment tensor Is a double couple. If It la not, tha moment tensor may be decom- 
posed either Into a double couple plus compensated linear dipole (Knopoff and Ran- 
dall, 1970) or two orthogonal double couples (GHbert, 1981). Generally, one double 
couple wOI.be larger than the other. Table A2 presents the major and minor double 
couples we obtained for our events. The notation conventions are the same as those 
In Kanamori and Given (1981), who follow Jarosch and Aboodl (1970). 

General moment tensor solutions often deviate significantly from double couples. 
Some investigators (Dziewonskl at aL, 1981; Dzlewonskl and Woodhouse, 1982; 
Giardlnl, 1982, 1983) have attached significance to these deviations. Dziewonskl et 
el. (1981) ’suggest a possible regional variation in the deviation from double couple. 
Giardlnl (ig82,1983) reports a possible depth dependence In the deviation, with high 
deviations occurring at intermediate depths In deeply penetrating subduction zones. 
The deviation from double couple Is a difficult parameter to interpret. There Is at 
present no way to evaluate the significance of variations in this parameter, end one 
cannot be certain that the data available can effectively constrain It (Kanamori and 
Given, 1982). 

The deuble couple is still a very useful model of the earthquake source, and one 
that has never been convincingly ruled out by the observations. Kanamori and Given 
(1981) present a technique for Inverting IDA data when the source is constrained to 
be a double couple. We have performed such an Inversion for each of our events and 
the results.are shown in Table A3. 

In general wc see that tho inversion of IDA data can be carried out qK*e suc- 
cessfully on deep and Intermediate events, despite the lower signal levels of surface 
waves as compared to shallow events. There Is, however, one difficulty which arises 
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at intermediate depths: tha ax citation functions of Rayleigh wavaa (saa Kanamori 
and Stewart, 1076) go through zeroes, as shown In Fig. A2. Kanamori and Given's 
matrix aquation (aquation (7) In tha 1081 papar) la 

AM * V 

whara M Is tha paramatar vactor, V is a vector of spactral points from tha data, and 
A is a matrix whosa antriss dapand on tl a ax citation functions. If any of tha exdta- 
tion functions vanish, A bacomas singular. When P R W vantshaa, tha aiamants M^ 
and ara Indatarminata. Whan Sg^ vanishes (M n +M m ) Is Indatarminata. 

Tha function Qg^ x) does not vanish axcapt at 2 aro dapth, so It Is not of any concam 
axcapt for vary shallow avants. Fig. A2(a) plots varaus dapth for pariods from 
1 00 - 274 saconds; this raprasants tha practically usabla ranga of pariods for thasa 
Inversions. Wo saa that tha zero occurs soma whara between 80 and 160 km depth. 
From fig. A2(b), wa saa that Pg^ goes through a zero somewhere between 130 and 
200 km depth. Thus IDA Inversions ara subject to numerical difficulties over a size- 
able portion of tha Intermediate dapth range. To soma extant, wa can avoid *he 
problem by choosing an appropriate period. This technique, however, can push tha 
zero away from the dapth of Interest by only a limited margin, and complications ara 
caused by tha fact that tha dapth is not known exactly. Reported dapths in tha 
NOAA catalog probably have error bars of at least ♦/■ 25 km. Tha hypocentral dapth 
Is not varied as a free parameter In tha Kanamort-GIven schema. 

One way to Improve solutions for events In the problematical depth range is to 
use first motion data In conjunction with the surface wave data. Kanamori (1982, 
personal communication) has developed a technique to do this. When 
+ =0, the standard expression for P-wave displacement In a 
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homogeneous whole-space can be written 

v(i*,jp) = jV av sln s t fc sin2p-(A/ vv -IV HI )(^ein 8 i k co82{p) 

HU m +li w )( j<1 -3cos 8 i fk ))-^f v> 8in2i lk slnp-tf as sln2i fc co8p 

where i* is takeoff angle and p is azimuth. By estimating u from amplitude data one 
can, in principle, determine the moment tensor. Because of the extreme scatter typi- 
cally observed In P wave amplitude data, however, the technique uses a very crude 
amplitude measure, with +1 for a clear compraeeional arrival, 0 for an ambiguous 
arrival, and -1 for a clear tensional one. After obtaining a guess at the moment ten- 
sor from this simple method, one then picks the maximum alement as a reference ele- 

ment Mr. In the case when vanishes, [IL&lMsL ^ constrained to have the 

same value in the surface wave inversion as is given by the body wave inversion. In 

M (M —M ) 

the case when PrM vanishes, both - 7 *- and — — — are constrained. This 

Ur Ur 

technique has been used on events 07/22/80 and 05/13/78 indicated In Table A1. 
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Tabla A1 
Uatof^Etrthcjij^ 


Data 

H 

M 

S 

Lat. 

Lon. 

Dapth 

E7V 

rnimm 

Raf. 

07 22 80 

07 

06 

23.0 

-20.3 

■ESI 

122 

6.1 

6.8 

1 

06 13 78 

07 

08 

46.2 

-14.6 

Waa 

160 

6.7 

- 

1 

03 16 78 

22 

04 

40.1 

26.4 

WmM 

263 

6.1 


1,3 

03 07 78 

02 

48 

38.4 

32.0 

■Eta 

442 

M 

1 

1,3,6 

04 28 81 

21 

14 

66.7 

-23.7 

180.0 

640 

l 


2 

08 18 78 

21 

31 

26.3 

41.8 

130.8 

688 

6.1 


1 

07 20 80 

21 

20 

04.0 

-17.8 

E Tjjjj 

681 

1 

VTS 

1 

04 24 78 

01 

46 

08.0 

-20.8 


688 




10 17 78 

06 

43 

03.0 

"8.6 


601 





IWftrtnMi 

(1) Qtordlnl (1943) 

(2) DzJawonakl and Wenzhou## (1083) 

(3) Oclawonakl at at. (IfJBI) 

(4) QNbart (1081) 

(8) Walok, personal communication (1 082) 

Tfco*® rafarancaa ara to Invaatlgatora hava atudlad thaaa a van t a Indapandantly from ua. in 
moat eaaaa, tha agraamant batwaan our aolutlona and thalra la axoalknt. Aaf. (1), tha aa yat unpubilahad 
atudy by Qtotvlnl (1083), and *af. (2) w^a a too tha aourcaa for tha momant tanaor In vara Iona of amaHar 
aarthquakaa from 1 077-1082, a urvayad In Saetlon 2. 
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Table AS 


Fault Plana Inversions 
(Standard Errors In Parentheses) 


Event 


6 

X 

9 

03 16 78 

0.063 

47.0 

178.2 

120.9 

(0.011) 

(119) 

(4.9) 

(4.8) 

03 07 78 

0.69 

24.1 

167.3 

48.6 

(0.06) 

(11.6) 

(9.9) 

(13.6) 

04 28 81 

0.18 

42.8 

-12.3 

76.2 

(0.019) 

(6.4) 

(3.7) 

(6.6) 

08 16 79 

0.14 

20.6 

136.3 

60.9 

(0.009) 

(6.7) 

(139) 

(16.6) 

07 20 80 

0.036 

80.4 

-46.8 

-33.7 

(0.006) 

(4.2) 

(11.3) 

(9.3) 

04 24 79 

0.024 

77.4 

268.2 

78.8 

(0.006) 

(6.6) 

(21.9) 

(10.9) 

10 17 79 

0.36 

66.6 

-1 34.0 


(0.03) 

(3.4) 

(4.6) 
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Figure A1. Focal machanlsma for avanta llatad In Tabla A1. For 06/13/78 and 
07/22/80, tha major double couple in the combined surface-wave/flrst- 
motlon Inversion la plotted. For the reat of the events, results of the fault 
plane Inversion are plotted. Circles on nodal planes indicate slip vectors. 
Filled squares are compression axes, and open ones are tension axes. 
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Flgure A2. Rayleigh wave excitation functions (a) and (b) as a function of 
depth. Number labelling each of the curves indicates the period In 
seconds. 
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Appendix B: A Brief Outline of Dlroetional Statistics 


When one's data are directions in space, one must use special methods to 
analyse them statistically. In two dimensions, for Instance, a unit vector with azimuth 
1* and one with azimuth 360° yield a vector with azimuth 180° if one simply aver- 
ages azimuths. This is clearly an Incorrect mean direction. The literature concerning 
directional data is extensive, and includes two full length monographs, one by Mardla 
(1972), and another by Batschelet (1981). Both are good books, but the second is 
particularly clear for the beginner. 

One simple and Intuitively pleasing way to average unit directions is to obtain 
the resultant, and take its direction as the mean direction (Watson, 1966). in the 
example above, this would give 0° as the mean direction. The resultant can also pro- 
vide an effective measure of dispersion. The closer a sample of vectors is grouped 
about a mean direction, the larger will be the resultant. These are measures one 
might logically choose for data whose exact distribution function were unknown, but 
which seemed to be grouped about a mean direction. 

These simple measures are also the ones applied rigorously in Fisherian statis- 
tics (Fisher, 1953), which have been used in paleomagnetism (e.g., Watson and Irv- 
ing, 1957; Watson, 1956a). The Fisher distribution, also known as the Von Mises dis- 
tribution on a sphere, is essentially a spherical analogue of the two dimensional Gaus- 
sian distribution. The density function is given by 
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6 and <p are, respectively, the polar and azimuthal angles, and k Is the Fisher preci- 
sion parameter. The mean Is, as we have said, of the same direction as the resultant, 
and the circular variance Is given by 


o*= 


N-R 

N 


-e 


k 


where N is the number of samples and R is the magnitude of the resultant vector; k is 
the maximum likelihood estimate of k. The confidence cone for the mean direction 
can be calculated, to significance (1 -a), from 


a= 


1 


(N-Rc) n 


-1 




+ ' W C 1) MilpL w-x - «)»-■+ 

2 R-Rc +4 


New terms are taken as the discontinuities R*N-2,N-4,N-6 ... are passed. Note that 
c =cos 6, where Q is the half -angle of the confidence cone. The Fisherian distribution 
Is circularly symmetric about the mean. It is also antipodally asymmetric, or unipolar. 
That is, a direction is distinct from its antipode. When dealing with data, such as 
earthquake compression or tension axes, where a direction and Its antipode are 
equivalent, one must take care to project those directions which would cluster about 
the antipode of the mean back to their antipodes. 

A more general distribution than that of Fisher is the Bingham distribution (Bing- 
ham, 19B4; 1974), which is now also finding use in paleomagnetism (Onstott, 1980). 
This is not a circularly symmetric distribution, so that elongate probability patches 
about a mean are allowed, it Is also antipodally symmetric. The best estimate of 
the mean direction Is not necessarily the direction of the resultant, but follows from a 
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moment of Inertia analysis (Mardla, 1072; Onstott, 1080). Bingham's density func- 
tion has the form 




1 k,(t »!)»♦*,(£ «,)») 

4ird(fc„fc,) 


where 


* n 0 0 

may bo evaluated asymptotically (Bingham, 1064) or numerically (Onstott, 1080). 4, 
and fc| are Bingham's "concentration parameters". As Onstott points out, the squar- 
ing of cos 6 reflects the antipodal symmetry of the distribution. Bingham (1064) was 
able to write the likelihood function of his distribution In terms of the moment-of- 
Inertla matrix. The moment of inertia, about a fixed axis U(x, y,a), of N points on a 
unit sphere, each point of unit mass and direction (i^m^n*), can be written (Mar- 
dla.1072) 


M = U r BU 


where B-I -T, and 


E 1 * 8 E*» rrti 

E**"* E^H* 1 ! E *< 8 


The distribution of the N points Is descrlbable by the eigenvalues of the matrix T. 
As discussed by Mardla and Onstott, Bingham's distribution function may reduce to 
the uniform distribution, or variously describe symmetric and asymmetric girdle distri- 
butions, or distributions about a mean direction (maximum eigenvector). When the 
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eigenvalues r< of T are distinct, with r 1 <T a «r 3 , the concentration parameters are 
such that ki<kz< 0 , and we have an elongate distribution about the maximum eigen- 
vector. When fci-»fc e , a circularly symmetric disributlon Is approached, and the max- 
imum eigenvector approaches the resultant vector. That is, we approach a situation 
similar to the Fisherian one. The variance in Bingham statistics is 

2 1 

where -tj). is the semiaxis of the standard deviation ellipse 

around the eigenvector in the direction of the 3 th eigenvector. Thus we are 
interested in o at and 032 , the semiaxes of the ellipse around the mean direction. The 
semi axes of the confidence region around the mean, to significance (1 -a), are given 
approximately by 



where x*,_ a ( 2 ) denotes the x 2 distribution with 2 degrees of freedom at percentage 
point ( 1 -a). 

Apart from performing statistical analyses based on a Fisher, Bingham, or other 
distribution, one can use the resultant vector of f t data in a test for randomness. 
By "randomness" here we really mean "uniformity" or "isotrcry", where no preferred 
direction exists. The null hypothesis is, then, that the parent population is uniformly 
distributed. One can compute the probability P (Watson, 1956b) that a uniform dis- 
tribution of vectors will yield a given resultant /? 0 . Hence, given a resultant exceed- 
ing Rq in magnitude, one can reject the hypothesis of uniformity, with confidence P. 
Good explanations and tables for performing the test are given in Stephens (1964; 



-233- 


1 960a, b). Variants of this test, as well as other tests, are discussed In Batschelet 
(1981). 



Appendix C: Plots of Seismlelty Versus Depth for the 
World's Subduetlon Zones 


In this appendix wa present histograms of the total niober of earthquakes with 
one-second body wave magnitude 4, for the time period 1004-1980, versus 
depth in the world's subduction zones. The data sources are the NOAA and PDE cata- 
logs. We note that we have regionalised such that the isolated areas of very deep 
seismicity in South America and the New Hebrides do not appear on the plots. 
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Flgure Cl. 


The world's subduction zones for which seismicity Is plotted In fig. C2 



r 


- 238 - 


ORIGINAL PAG* IS 
OF POOR QUALITY 



Figure C2. Plots of log^tf versus depth for the regions in fig. Cl. AT Is the total 
number of events (1964*1 980) with m^>4. 
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Appendix Dt Spatial Stress Plots and Velocity Plots not 
Included In the T»xt 


We present here spatial stress plots and flow fields not included In the text, as 
discussed in Section 3.1.2. Figures are labelled according to the following code: 
Stress plots are labelled DnnS, and velocity plots are labelled OnnV. "0" simply 
Indicates Appendix D, S and V respectively indicate that the figure is a stress plot or 
a velocity plot, and nn Is the numjor of the figure in if le taxt to which the particu- 
lar figure in this appendix corresponds. 
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